IMPLEMENTATION  OF  BRANCH-POINT-TOLERANT  WAVEFRONT 
RECONSTRUCTOR  FOR  STRONG  TURBULENCE  COMPENSATION 

THESIS 


Michael  J.  Steinbock,  Civ,  USAF 
AFIT/GE/EN  G/ 1 2-45 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Patterson  Air  Force  Base,  Ohio 


DISTRIBUTION  STATEMENT  A: 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  United  States  Air  Force,  the  Department  of  Defense,  or  the  United 
States  Government. 

This  material  is  declared  a  work  of  the  U.S.  Government  and  is  not  subject  to  copyright 
protection  in  the  United  States. 


AFIT/GE/ENG/12-45 


IMPLEMENTATION  OF  BRANCH-POINT-TOLERANT  WAVEFRONT 
RECONSTRUCTOR  FOR  STRONG  TURBULENCE  COMPENSATION 


THESIS 


Presented  to  the  Faculty 

Department  of  Electrical  and  Computer  Engineering 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
in  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science  in  Electrical  Engineering 


Michael  J.  Steinbock,  B.S.E.E.,B.S.Co.E. 
Civ,  USAF 

June  2012 


DISTRIBUTION  STATEMENT  A: 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AFIT/GE/ENG/ 12-45 


IMPLEMENTATION  OF  BRANCH-POINT-TOLERANT  WAVEFRONT 
RECONSTRUCTOR  FOR  STRONG  TURBULENCE  COMPENSATION 


Michael  J.  Steinbock,  B.S.E.E.,B.S.Co.E. 
Civ,  USAF 


^ .  a,,!- 


J^SQhjD.  Schmidt,  Maj,  US AFR  (Committee  Member) 

L 


(Committee  Member) 


Si  201-9- 


Date 

\ 

ifin 

Date 

SJu 

Da1^ J 

AFIT/GE/ENG/12-45 


Abstract 


Branch  points  arise  in  optical  transmissions  due  to  strong  atmospheric  turbulence,  long 
propagation  paths,  or  a  combination  of  both.  Unfortunately,  these  conditions  are  very  often 
present  in  desired  operational  scenarios  for  laser  weapon  systems,  optical  communication, 
and  covert  imaging,  which  suffer  greatly  when  traditional  adaptive  optics  systems  either 
cannot  sense  branch  points  or  implement  non-optimal  methods  for  sensing  and  correcting 
branch  points.  Previous  research  by  Pellizzari  presented  a  thorough  analysis  of  various 
novel  branch  point  tolerant  reconstructors  in  the  absence  of  noise.  In  this  research  a  realistic 
model  of  the  Air  Force  Institute  of  Technology’s  adaptive  optics  system  is  developed  and 
used  for  the  first  realistic  tests  of  these  branch  point  tolerant  reconstructors.  Utilizing 
a  self-referencing  interferometer  as  the  high-order  wavefront  sensor-this  type  of  sensor 
being  theoretically  immune  to  scintillation-this  effort  has  extended  previous  theoretical 
work  by  adding  realistic  noise  effects  to  the  SRI’s  measurements  before  reconstructing  the 
wavefronts  and  applying  the  control  law.  Not  only  is  the  adaptive  optics  correction  shown 
to  outperform  the  traditional  techniques  by  as  much  as  126%,  but  several  new  theories 
and  refinements  to  existing  theories  were  discovered  along  the  way.  This  study  provides 
a  foundation  to  guide  hardware  implementation  in  the  future,  where  noise  effects  will  be 
present. 
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IMPLEMENTATION  OF  BRANCH-POINT-TOLERANT  WAVEFRONT 
RECONSTRUCTOR  FOR  STRONG  TURBULENCE  COMPENSATION 

I.  Introduction 

Within  the  Department  of  Defense  (DoD),  directed  energy  (DE)  capabilities 
have  been  identified  as  one  of  “the  key  ‘game-changing’  technology-enabled 
capabilities  that  enters  service  during  [the  next  two  decades]”  [44].  DE  spans  a  wide 
range  of  potential  applications  yet  has  many  technological  challenges  still  to  be  overcome. 
Applications  of  particular  interest  to  the  DoD  are  space  situational  awareness,  laser 
communication,  and  laser  weapon  systems  [44].  Currently,  a  major  obstacle  in  each 
of  these  applications  is  accounting  for  the  optical  aberrations  induced  by  the  random 
nature  of  Earth’s  atmosphere.  For  a  point  of  reference,  these  same  atmospheric  effects 
are  what  prompted  NASA  to  spend  billions  of  dollars  to  design  and  maintain  the  Hubble 
Space  Telescope.  To  deal  with  these  atmospheric  aberrations,  there  are  only  a  handful  of 
possibilities.  Adaptive  optics  (AO)  utiliy.es  a  wavefront  sensor  (WFS)  and  reconstructor  to 
estimate  the  optical  distortions  on  some  reference  beam  traveling  through  the  atmospheric 
path  of  interest.  Subsequently,  a  wavefront  corrector,  such  as  a  deformable  mirror 
(DM),  removes  the  estimated  distortions.  While  there  are  other  methods  for  dealing  with 
Earth’s  random  effects  on  light,  AO  is  unique  in  providing  real-time,  improved  imaging 
performance,  which  is  essential  for  the  above  applications  of  interest. 

Current  AO  techniques  only  provide  good  correction  capability  for  limited  operational 
regimes.  As  the  optical  propagation  paths  lengthen  and/or  the  atmospheric  turbulence 
becomes  stronger,  the  propagated  light  begins  to  constructively  and  destructively  interfere 
with  itself  [5].  In  these  circumstances,  traditional  wavefront  reconstruction  techniques 
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used  to  drive  the  optical  correction  process  are  unable  to  completely  estimate  the  optical 
aberrations.  These  situations  are  of  definite  operational  interest  to  the  DoD,  but  most 
current  AO  techniques  lose  effectiveness  in  these  more  challenging  situations  [31,39,42]. 
The  research  presented  here  is  aimed  at  extending  the  operational  regime  of  AO  to 
offer  improved  correction  capability  when  optical  interference  begins  to  severely  degrade 
performance  of  traditional  AO  techniques. 

1.1  Problem  Statement 

Branch  points  arise  in  optical  transmissions  due  to  strong  atmospheric  turbulence, 
long  propagation  paths,  or  a  combination  of  both.  Unfortunately,  these  conditions  are 
very  often  present  in  desired  operational  scenarios  for  laser  weapon  systems,  optical 
communication,  and  covert  imaging,  which  suffer  greatly  when  traditional  adaptive 
optics  systems  either  cannot  sense  branch  points  or  implement  non-optimal  methods  for 
estimating  and  correcting  branch  points.  The  overall  goal  for  this  research  is  to  transition 
previous  theoretical  work  in  optimized  wavefront  reconstruction  for  strong  turbulence 
AO  to  an  initial  implementation  in  hardware.  To  bridge  the  transition,  a  high-fidelity 
computer  model  is  developed  and  used  in  testing  the  new  algorithms  prior  to  hardware 
implementation.  The  work  in  simulation  is  shown  to  be  invaluable  in  further  investigating 
the  behavior  of  the  new  wavefront  reconstruction  algorithms  with  many  unanticipated 
lessons  and  discoveries. 

1.2  Key  Results 

This  research  primarily  investigated  the  behavior  and  properties  of  post-processing 
congruence  operation  (PCO)  based  wavefront  reconstructors.  The  following  list  summa¬ 
rizes  the  contributions  made  in  this  research  towards  bringing  the  PCO-based  reconstructor 
to  an  implementable  state. 
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•  Demonstrated  greatly  improved  strong  turbulence  AO  performance  in  terms  of  both 
mean  Strehl  and  Strehl  variance.  The  AO  performance  gains  using  a  PCO  wavefront 
reconstructor  over  the  traditional  least  squares  reconstructor  in  terms  of  mean  field 
estimated  Strehl  ratio  and  Strehl  ratio  normalized  standard  deviation  are  summarized 
in  Table  1.1.  Mean  Strehl  is  a  metric  of  absolute  AO  correction  performance  (higher 
is  better),  and  the  Strehl  standard  deviation  is  a  metric  of  correction  stability  (lower 
is  better).  These  results  were  presented  at  the  2012  IEEE  Aerospace  conference. 


Table  1.1:  Aggregate  AO  performance  gains  of  all  PCO-based  wavefront  reconstructors 
over  the  traditional  least  squares  reconstructor. 


Mean  Strehl 

Normalized  Strehl  Standard  Deviation 

0.04 

0.2% 

-2.4% 

0.5 

19.1% 

-50.6% 

1 

126.4% 

-57.6% 

•  Identified  cause  of  fades  from  PCO-based  AO  and  demonstrated  a  potential  solution. 
When  in  closed-loop  operation  in  strong  turbulence,  PCO  algorithms  can  cause 
short-term  fades  as  hopt  varies.  These  drops  in  Strehl  are  not  caused  by  the  higher- 
order  control  law  as  much  as  different  tilt  components  associated  with  each  PCO 
rotational  realization.  The  tilt  is  removed  from  the  higher-order  corrections  and 
instantly  dumped  on  the  tracking  system;  causing  the  fades.  One  solution  was  found 
by  allowing  tilt  in  the  higher-order  corrections  with  conservative  controller  gains. 
This  prevented  the  fades  with  a  slight  tradeoff  in  mean  Strehl  ratio.  Unfortunately,  a 
realistic  system  with  slight  misalignments  would  require  a  zero-DC-gain  control  law 
to  implement  a  similar  solution.  Additionally,  this  puts  higher  stroke  requirements 
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on  the  DM  when  tilt  is  being  included.  With  this  new  understanding  future  work  in 
control  design  could  help  mitigate  these  fades  based  on  the  application  requirements. 

•  Demonstrated  and  developed  theory  on  improved  PCO  performance  in  weak 
turbulence,  low  SNR  regimes.  In  weak  turbulence,  where  no  branch  points  are 
present,  a  consistent  performance  improvement  of  the  PCO  reconstructors  over 
traditional  least  squares  was  found.  Although  the  difference  is  small,  this  was 
completely  unexpected  based  on  previous  concerns  of  noise  corrupting  the  PCO 
process.  A  new  theory  proposed  by  the  author  shows  in  high-noise  conditions, 
least  squares  reconstructors  can  introduce  phase  anomalies  to  the  AO  corrections. 
However,  the  PCO  process  removes  these  phase  discrepancies  introduced  by  the  least 
squares  reconstruction.  These  findings  will  be  presented  at  the  2012  SPIE  Optics  and 
Photonics  conference. 

•  Developed  a  set  of  theories  regarding  phase  anomalies  in  high  wavefront  gradient 
conditions.  When  localized  wavefront  phases  being  reconstructed  from  an  SRI  are 
greater  than  n  radians  between  adjacent  subapertures,  the  least  squares  reconstructor 
will  alias,  leading  to  a  phase  anomaly  in  the  reconstructed  field.  This  can  be  mitigated 
using  a  higher  density  of  SRI  subapertures,  at  the  cost  of  SNR.  Alternatively,  initial 
results  show  the  PCO  corrects  the  phase  anomalies  introduced  by  the  least  squares 
reconstruction. 

•  Formulated  an  optimized  data  flow  for  an  AO  system  using  PCO.  Started  by 
diagnosing  coupling  between  the  high-order  and  tracking  control  systems  when  using 
a  PCO-based  wavefront  reconstructor.  Since  PCO  reconstruction  is  not  a  single 
matrix  operation,  decoupling  tilt  from  the  higher-order  commands  required  a  new 
data  flow  process.  The  resulting  process,  illustrated  in  Fig.  3.4,  minimizes  matrix 
multiplies  while  maximizing  AO  effectiveness. 
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•  Conducted  first  beam  projection  simulations  using  PCO  wavefront  reconstruction 
techniques.  Conducted  at  the  end  of  this  research,  this  establishes  a  foundation 
for  future  efforts  to  investigate  these  novel  reconstruction  methods  with  respect 
to  beam  projection  applications.  Although  the  initial  simulation  results  identified 
needed  changes,  initial  observations  indicate  promise  over  traditional  least  squares 
reconstruction.  The  identified  flaws  have  since  been  fixed  and  follow-up  simulation 
results  will  be  presented  at  the  2012  DEPS  Beam  Control  Conference. 

•  Developed  an  extensible  AFIT  AO  system  model  for  future  research.  The  model 
is  configurable  to  varying  levels  of  realism  and  AO  capabilities.  Allows  isolation 
of  the  research  problem  at  hand  without  the  difficulties  of  hardware  work  and  the 
assumptions  of  analytical  work.  This  includes  the  first  modeling  of  realistic  fiber 
coupling  to  our  specific  SRI  optical  setup  -  crucial  to  realistic  noise  modeling. 

1.3  Thesis  Overview 

Chapter  II  provides  an  overview  of  traditional  AO  and  the  optical  effects  as  a  result 
of  the  randomized  atmosphere.  This  is  followed  by  an  introduction  to  branch  points  and 
the  previous  theoretical  work  in  developing  an  optimized  branch-point  tolerant  wavefront 
reconstructor.  Chapter  III  describes  a  high-fidelity  computer  simulation  developed  to 
model  AFIT’s  AO  system.  This  model  exposed  new  insights,  crucial  to  a  successful 
implementation  of  the  new  reconstruction  techniques.  These  insights,  discussed  in  Ch.  IV, 
offer  explanations  to  a  number  of  previously  posed  questions,  and  a  new  extension  of  the 
theory  is  presented  for  beam  projection  applications.  Finally,  Ch.  V  presents  a  summary  of 
the  challenges  overcome,  key  contributions,  and  suggested  future  work  from  this  research 
effort. 


5 


II.  Background  and  Related  Research 


This  chapter  introduces  the  current  state  of  AO  and  previous  theoretical  methods 
for  extending  AO  operability  into  regimes  with  scintillation  present.  An  overview  of 
a  typical  AO  system  is  given,  along  with  two  WFS  types  -  each  with  unique  strengths 
and  weaknesses.  Following  this,  the  detrimental  effects  of  scintillation  are  explained 
with  respect  to  AO.  This  chapter  concludes  by  discussing  previous  research  in  specialized 
wavefront  reconstruction  techniques  that  attempt  to  improve  upon  traditional  reconstructors 
when  branch  points  are  present  in  the  received  beacon  wavefronts. 

2.1  Traditional  Adaptive  Optics 

As  light  propagates  through  Earths  atmosphere,  stochastic  variations  in  temperature 
and  pressure  cause  the  light  to  both  diffract  and  refract  [5].  This  causes  beam  wander,  beam 
spreading,  and  phase  fluctuations.  A  simple  example  of  this  can  be  observed  in  the  apparent 
twinkle  of  stars  at  night.  The  Hubble  Space  Telescope  was  built  to  operate  beyond  the 
random  atmosphere,  where  the  stars  appear  as  sharp  objects.  For  many  applications,  these 
atmospheric  effects  severely  limit  performance.  For  instance,  beam  projection  systems 
rely  on  efficient  delivery  of  light  to  a  small  target  area,  and  optical  communication  must 
maintain  a  reliable  optical  link.  One  remedy  is  to  estimate  the  atmospheric  distortions 
and  correct  for  them  in  real  time.  This  approach  is  called  AO,  utilizing  a  WFS  and 
reconstructor  to  estimate  the  optical  distortions  on  some  reference  beam  followed  by  a 
wavefront  corrector  to  mitigate  the  estimated  distortions. 

2.1.1  Atmospheric  Turbulence. 

Within  the  atmosphere,  light  is  distorted  as  it  encounters  spatially  and  temporally 
randomized  variations  in  the  air’s  index  of  refraction.  The  statistical  behavior  of  these 
random  variations  were  first  described  in  the  1940’s  by  Kolmogorov,  when  he  used  physical 
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Figure  2.1:  Turbulent  eddies  shown  of  size  L0  down  to  Iq. 


insights  to  describe  the  velocity  of  air  within  the  atmosphere  [5].  Energy  cascade  theory 
says  that  kinetic  energy  injected  to  a  viscous  fluid  creates  large  sized  turbulent  eddies. 
These  large  eddies  break  down  into  progressively  smaller  eddies  before  the  energy  is 
eventually  transferred  into  heat.  The  dimensions  of  these  eddies  are  described  by  the 
inertial  subrange,  where  eddie  sizes  fall  between  the  inner  scale,  /0,  and  outer  scale,  L0. 
Within  the  inertial  subrange,  the  turbulence  can  be  considered  statistically  homogenous 
and  isotropic  [5].  Figure  2.1  provides  a  visual  depiction  of  the  energy  cascade  theory  as 
kinetic  energy  is  converted  down  into  thermal  energy. 
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With  this  physical  foundation,  Kolmogorov  was  able  to  formulate  a  one  dimensional 
structure  function  for  wind  velocities  within  the  inertial  subrange  [5].  From  this,  one 
can  derive  three  dimensional  power  spectrum  of  the  velocity  fluctuations  within  the 
atmosphere  [17]: 

Okr(k)  =  0.066e2/V11/3,  (2.1) 

where  k  is  the  scalar  spatial  frequency  (in  units  of  rad/m)  and  e  is  the  average  energy 
dissipation  rate  of  the  atmosphere.  By  defining  C\  as  the  velocity  structure  constant, 

C2V  =  2e2/3,  (2.2) 


Eq.  (2.1)  can  be  simplified  to 

Or*(k)  =  0.033C2/T11/3.  (2.3) 

The  velocity  structure  constant,  with  units  of  m4/3/s2,  provides  a  quantitative  evaluation  of 
the  energy  within  the  atmospheric  turbulence  and  is  actually  based  in  the  formulation  of 
the  velocity  structure  function  [5]. 

From  Kolmogorov’s  work,  Obukhov  [28]  and  separately  Corrsin  [8]  were  able 
to  extended  his  statistical  model  of  atmospheric  velocity  variations,  first  to  potential 
temperature,  and  later  to  index  of  refraction  variations.  This  important  final  leap  directly 
relates  the  phase  aberrations  accumulated  by  light  to  the  randomized  atmosphere.  The 
resulting  three  dimensional  refractive  index  power  spectrum,  known  as  the  Kolmogorov 
spectrum,  is  given  by 

$>„(*)  =  0.033 Cyn/3,  (2.4) 

where  n  is  the  index  of  refraction  and  C2  is  the  index  of  refraction  structure  constant.  C2  is 
a  physical  quantity  that  can  be  measured  for  a  given  location,  time  of  day,  and  altitude. 
Several  models  of  C2  exist  that  can  be  used  for  simulations  of  a  desired  scenario  [5]. 
Figure  2.2  shows  the  commonly  used  Hufnagel- Valley  model,  along  with  many  of  the 
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Figure  2.2:  The  Hufnagel- Valley  C\  profile  model.  Figure  from  [23]. 


physical  atmospheric  layers  that  drive  its  behavior.  For  this  research,  a  horizontal  path 
was  assumed  -  implying  a  constant  C\  along  the  propagation  path. 

It  is  important  to  note  the  power  spectral  models  for  velocity,  temperature,  and 
index  of  refraction  variations  all  follow  similar  behavior  -  easily  seen  by  comparing 
Eqs.  (2.3)  and  (2.4).  However,  as  presented  thus  far,  effects  outside  the  inertial  subrange 
have  been  ignored.  The  Kolmogorov  spectrum  actually  assumes  an  infinitely  large  outer- 
scale  and  infinitely  small  inner-scale.  This  is  an  assumption  made  for  mathematical 
simplicity,  but  does  not  match  the  true  atmospheric  fluidics.  Hill  developed  a  numerical 
model  based  on  a  more  detailed  hydrodynamic  analysis  [20],  later  adapted  to  analytic  form 
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by  Andrews  [4],  given  by 


<&„(*)  =  0.033C" 

0  <  K  <  OO, 


1  +  1.802 


2541  — 


7/6 


exp  k2/kj J 


N  -  3.3 //o, 


K,l  I  (k2  + 

/Co  =  2n/L0. 


11/6  ’ 


(2.5) 

(2.6) 


This  is  sometimes  called  the  modified  atmospheric  spectrum,  where  the  math  to  the  right 
of  C\  has  been  added  to  Eq.  (2.4)  to  match  the  results  found  by  Hill.  Most  notably,  this 
spectrum  model  includes  a  spectral  bump  at  wavenumbers  near  the  inner-scale  size.  As 
the  inner-scale  size  decreases,  light  begins  to  diffract  at  the  smaller  eddies  [5].  With  the 
Hill  spectrum’s  increased  content  of  small  eddie  sizes,  scintillation  effects  are  seen  earlier 
than  a  similar  propagation  using  a  Kolmogorov  spectrum.  For  this  research,  the  modified 
spectrum  was  used  in  all  numerically  propagated  optical  fields. 

Mathematically,  Maxwell’s  Equations  govern  all  electromagnetic  interactions.  In 
optics,  these  equations  are  typically  simplified  by  making  three  key  assumptions: 

1.  Negligible  backscattering  and  polarization  effects. 

2.  The  Markov  approximation.  Assumes  the  atmospheric  refractive  index  variations  are 
delta  correlated  in  the  direction  of  propagation. 

3.  The  paraxial  approximation.  Assumes  all  light  rays  of  interest  travel  at  angles  very 
near  the  optical  axis. 

From  these  simplifications,  several  methods  exist  for  applying  Maxwell’s  equations  to 
propagating  an  initial  field  through  a  random  medium.  For  regimes  of  weak  turbulence, 
one  such  method  is  the  Rytov  approximation  [5]: 

U( r,  L)  =  U0(r,  z  =  L)  exp  [T(r,  L)] ,  (2.7) 

where  the  initial  field  is  propagated  a  distance  L,  t/0(r,  L)  is  the  field  propagated  through 
vacuum,  and  it  is  modified  by  the  complex  phase  perturbation,  T^r,  L).  Further,  the 
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perturbation  can  be  split  into  amplitude  and  phase  components  [38] 


'F  =  X  +  i(P- 


(2.8) 


where  p  is  the  distance  from  the  optical  axis,  k  =  2n/A  is  the  optical  wavenumber,  J0 
is  a  Bessel  function  of  the  first  kind,  order  zero,  and  0„(/c)  is  the  atmospheric  index 
of  refraction  spectrum  model  being  used.  For  analytic  tractability,  the  Kolmogorov 
atmospheric  spectrum  is  assumed  for  the  following  definitions. 

By  substituting  the  Kolmogorov  spectrum  into  Eq.  (2.1 1),  further  substituting  that  into 
Eq.  (2.10),  and  solving  for  the  e~x  point  one  can  find  the  spatial  coherence  radius  of  the 
field,  p0.  However,  the  more  commonly  used  term,  r0  =  2.1po,  called  the  seeing  parameter, 
can  be  shown  to  be  [5] 

rL  /7\5/ 3  r3/5 

ro,sPh  =  OA13k2Jo  dz  ■  (2.12) 
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This  is  also  referred  to  as  Fried’s  parameter  after  being  introduced  by  D.  L.  Fried  [9]. 
With  respect  to  AO,  r0  is  a  critical  parameter  for  designing  a  system.  Apertures  larger  in 
diameter  than  r0  require  AO  to  see  additional  gains  in  imaging  resolution.  Further,  the 
spatial  density  of  wavefront  estimation  and  correction  by  an  AO  system  are  often  related  to 
tq.  To  be  effective,  roughly  one  WFS  subaperture  and  DM  actuator  spacing  should  span  r0. 

To  describe  the  temporal  coherence  of  the  field  aberrations,  Greenwood  analyzed  an 
ideal  AO  system  with  a  single  pole  control  filter.  He  calculated  the  characteristic  frequency 
of  the  system  with  respect  to  the  ratio  of  aperture  diameter  to  r0  [19]  and  derived  what  is 
now  called  the  Greenwood  frequency: 

[  rL  i3/5 

fg  =  0.254fc6/5  C2n(z)v5l\z)dz  ,  (2.13) 

Jo 

where  v(z)  is  the  perpendicular  wind  velocity  along  the  propagation  path. 

2.1.2  Wavefront  Sensing. 

The  first  step  to  mitigate  the  atmospheric  effects  on  light  is  to  be  able  to  estimate  the 
distortions  introduced  for  a  given  propagation  path  and  instant  in  time.  This  process  is 
called  wavefront  sensing  and  requires  some  form  of  reference  beacon  transmitting  light 
through  the  turbulence  column  of  interest.  With  longer  electromagnetic  wavelengths,  such 
as  in  radar,  the  fields  of  interest  can  be  directly  measured  and  compared.  Unfortunately,  in 
optics  the  short  wavelengths  and  corresponding  high  frequencies  are  beyond  the  detection 
range  of  current  electronics.  In  optics,  wavefront  estimation  requires  inferring  information 
about  the  received  beacon  field  from  other  physical  processes  [41].  The  following 
subsections  give  an  overview  of  three  such  methods. 

2.1.2. 1  Tilt  Tracker. 

One  of  the  most  simplistic  wavefront  sensors  is  a  tilt  sensor.  Shown  in  Fig.  2.3,  this 
sensor  is  meant  to  estimate  the  low-order  wavefront  aberrations  of  x  and  y  tilt.  By  focusing 
incoming  light  onto  a  quadcell  detector  or  sensor  array,  the  deflection  of  the  focused  light 
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can  be  calculated  by  the  centroid  equations 


Ax  = 


y  xi 

Zjx= 1 
yN  , 
Zjx=1  1x 


Ay  = 


■y. 


(2.14) 


yl\  , 

Zjy=l  ly 

where  x  and  y  denote  detector  column  and  row  positions,  x  and  y  are  the  mean  column 
and  row,  and  Ix  and  Iy  are  the  summed  detected  irradiance  values  for  each  column  and  row. 
Using  the  results  from  Eq.  (2.14),  the  incoming  wavefront’s  mean  tilt  can  be  determined  as 


ex  =  Ax/f,  9y  =  Ay/  f,  (2.15) 

where  /  is  the  focal  length  of  the  focusing  optic  used,  shown  in  Fig.  2.3. 

Tilt  comprises  approximately  87%  of  all  atmospheric  aberrations  [27],  and  as  such, 
using  a  dedicated  fast-steering  mirror  (FSM)  to  compensate  tilt  alone  offloads  a  large  stroke 
requirement  from  the  high-order  wavefront  corrector  [25].  Further,  the  tracker  can  be  used 
with  broadband  light,  has  no  moving  parts,  and  is  computationally  simple  to  operate.  These 
reasons  make  the  tracker  a  very  robust  wavefront  sensor;  however,  it  can  only  sense  the 
lowest  order  wavefront  aberrations. 


2. 1.2.2  Shack-Hartmann  Wavefront  Sensor. 

The  Shack-Hartmann  WFS  uses  a  tiled  array  of  trackers  across  the  aperture.  This 
divides  the  full  aperture  into  subapertures,  whose  incident  wavefronts  are  assumed  to 
be  roughly  planar.  With  this  architecture,  local  wavefront  tilts,  or  gradients,  are  sensed 
which  are  then  pieced  together  using  a  wavefront  reconstruction  algorithm  -  discussed  in 
Sec.  2.1.3.  Because  of  this  aperture  averaging  effect  with  respect  to  high  order  wavefront 
aberrations,  it  is  advisable  to  have  at  least  one  subaperture  per  r0  spacing  of  the  wavefront 
being  sensed  [25]. 

Figure  2.4  shows  a  Shack-Hartmann  design  using  square  lenslets  that  focus  down  on  to 
2x2  pixel  regions  on  a  larger  focal  plane  array.  This  is  a  typical  design,  but  other  variations 
exist  for  a  range  of  applications.  One  common  design  decision  is  to  include  a  band  of 
unsensed  pixels  between  subapertures  to  prevent  light  from  one  subaperture  ‘spilling’  into 
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Figure  2.3:  Principles  of  operation  for  a  tilt  sensor  where  6  is  the  incoming  wavefront  mean 
tilt  with  respect  to  the  sensor’s  optical  axis,  d  is  the  lens  diameter,  /  is  the  effective  focal 
length,  and  Ax  is  the  lateral  shift  of  the  focal  pattern  irradiance. 


an  adjacent  subaperture’s  sensed  pixels.  However,  this  requires  the  use  of  a  larger  and  more 
expensive  focal  plane  array  and  typically  increases  the  required  sensor  readout  time.  This 
wavefront  sensor  is  very  widely  used  as  it  shares  the  same  benefits  as  the  tracker  does. 
The  primary  downside  of  this  WFS  is  its  high  sensitivity  to  alignment  errors.  As  in  most 
WFS’s  there  is  a  tradeoff  to  be  made  between  dynamic  range  of  detectable  wavefront  tilts, 
measurement  accuracy,  signal  sensitivity,  and  speed. 

2. 1.2.3  Self-Referencing  Interferometer  Wavefront  Sensor. 

The  Self-Referencing  Interferometer  (SRI)  WFS  uses  interference  effects  to  directly 
sense  an  incoming  field,  as  opposed  to  most  WFS  designs  which  measure  some  form  of 
wavefront  gradients.  It  is  based  on  the  point  diffraction  interferometer  (PDI),  shown  in 
Fig.  2.5.  The  PDI  focuses  all  incoming  light  to  a  partially  transmissive  filter  in  the  focal 
plane.  The  on-axis  point  in  the  filter  is  completely  transmissive,  acting  as  a  spatial  filter. 
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Figure  2.4:  Diagram  of  the  Shack-Hartmann  WFS. 


When  the  light  is  recollimated,  an  interferogram  between  the  aberrated  input  wavefront  and 
spatially  filtered  reference  wavefront  can  be  seen  [24, 34] . 

Unfortunately,  the  single  interferogram  produced  by  a  PDI  does  not  give  a  direct 
means  of  forming  a  wavefront  estimate,  but  this  is  where  the  SRI  extends.  Here,  the 
incoming  light  is  split  into  two  separate  paths  before  being  recombined.  The  first  path 
is  responsible  for  spatially  filtering  the  incoming  wavefronts,  similar  to  the  pinhole  filter 
used  in  the  PDI.  However,  in  the  SRI  a  single-mode  fiber  is  used  instead  of  a  pinhole  filter. 
Differently  from  the  PDI,  this  filtered  reference  beam  is  further  split  up  into  separate  beams, 
and  each  beam  receives  a  different  phase  shift  -  explained  below.  The  second  path  applies 
no  filtering  or  altering  of  the  beam,  but  is  carefully  matched  in  length  to  the  reference 
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Figure  2.5:  Diagram  of  the  Point  Diffraction  Interferometer. 


legs  to  maintain  a  high  degree  of  temporal  coherence.  Finally,  each  phase  shifted  beam  is 
recombined  with  the  beacon  wavefront  to  form  an  interferogram  for  each.  The  irradiance 
pattern  of  each  interferogram  is  measured  on  an  imaging  sensor  for  use  in  calculating  the 
wavefront  phase.  Figure  2.6  shows  a  simple  layout  for  an  SRI. 

A  direct  measurement  of  the  input  wavefront  phase  can  be  calculated  using  the 
separate  interferograms  [11,34].  Assuming  a  perfectly  coherent  source,  the  m-th  sensed 
pixel  of  the  n- th  interferogram  follows 


h(- xm,ym)  =  I  (  \Ub(x,y)  +  Ur(x,y)e  'e"\~  dx dy,  (2.16) 

Am 

where  Ub  is  the  input  beacon  field,  Ur  is  the  reference  field,  and  Am  is  the  area  of  the 
corresponding  pixel.  If  we  assume  the  amplitudes  of  both  fields  are  equal  to  A/2,  Eq.  (2.16) 
can  be  shown  to  equal 


7n(xm,ym)  A  [l  +  COS i^(pbxm,ym  "b  @n) j  •  (2.1V) 

Here,  the  planar  nature  of  the  reference  wavefront  is  clear,  while  the  aberrated  input 
wavefront  remains  a  function  over  xm  and  ym.  There  are  multiple  methods  for  calculating 
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Figure  2.6:  The  Self-Referencing  Interferometer. 


the  input  wavefront  phase  using  different  numbers  of  interferograms  with  different 
reference  phase  shifts  applied  [11].  For  this  research,  a  four  bin  approach  is  used.  In 
this  configuration,  phase  shifts  of  0,  /r/2,  n,  and  3n/2  are  applied  respectively  to  each 
reference  beam,  each  interferogram  constituting  a  bin.  Following  Eq.  (2.17),  the  measured 
irradiances  of  each  interferogram  become 


(-^m>  ym) 

A  +  A  cos (cpb( 

ii 

5 

[l  +  COS  (ct>b(xm,ym)  +  f ) 

|  =  A  -  A  sin(06( 

ym)  ~  A  | 

[l  +  cos  {(pb{Xm,ym)  +  7t)] 

=  A  -  A  cos((pb{_ 

J-4  ^nv>ym)  A\_ 

1  +  COS  (fa(Xm,ym)  +  y) 

]  =  A  +  A  sin(0fc( 

(2.18) 


From  this  system  of  equations,  the  aberrated  input  wavefront  can  be  easily  calculated 
according  to  Eq.  (2.19)  and  (2.20): 


tan  ((pb(Xm,ym)) 


sin ((pb(xm,ym))  _  h(xm,ym)  ^2(xm,ym) 
cos (4>b(,xmiym))  h(xm,ym)  ~  h(xm,ym) 


(2.19) 
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(i PsRI  =  Tan  1 


(2.20) 


For  some  applications,  including  the  strong  turbulence  reconstruction  techniques  in 
Sec.  2.3,  it  is  important  to  note  that  the  complete  field  estimate  is  available  with  the 
relationship 

UsR^ih-hl  +  Uh-h),  (2.21) 

from  which  the  field  irradiance  can  be  calculated. 

Equation  (2.20)  reveals  a  remarkable  point  of  interest  for  real-time  AO:  using  an 
SRI,  the  wavefront  phase  can  be  directly  measured  by  merely  subtracting  the  measured 
interferograms,  one  division,  and  performing  an  inverse  tangent  operation  on  that  result. 
Further,  for  a  spatial  shifting  SRI  these  operations  can  be  performed  independently  for 
each  subaperture,  meaning  that  calculations  can  begin  immediately  as  pixel  data  from  each 
interferogram  is  read  in  from  the  sensor(s).  To  capitalize  on  this  advantage,  thought  must 
be  given  to  the  pixel  readout  order  of  the  sensor(s).  The  computational  time  required  by 
the  arctangent  operation  is  a  tiny  fraction  of  the  full  wavefront  reconstruction  required  for  a 
Shack-Hartmann.  However,  the  output  of  the  arctangent  is  wrapped  |-/r,  n\.  which  is  well 
suited  for  a  segmented  DM,  but  is  difficult  to  reproduce  over  a  continuous  face-sheet  DM.  In 
the  case  of  a  continuous  face-sheet  DM,  one  needs  to  implement  a  wavefront  reconstruction 
algorithm  in  the  control  system  [10, 14, 30].  Due  to  the  influence  functions  of  a  continuous 
face  DM,  sharp  wrapping  cuts  are  difficult  to  reproduce  and  may  even  cause  damage  to 
hardware.  Unfortunately,  the  reconstruction  process  eliminates  the  efficiency  gains  seen 
over  the  Shack-Hartmann  WFS. 

The  SRI  still  has  one  major  advantage  over  other  WFS  technologies  in  that  it  is  nearly 
immune  to  the  effects  of  scintillation  seen  in  strong  atmospheric  turbulence.  A  scintillated 
field  has  null  points  in  irradiance  over  the  phase  wavefronts  that  are  being  estimated  by  a 
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wavefront  sensor.  At  these  points  of  zero  irradiance  field  estimation  issues  arise  in  other 
wavefront  sensing  technologies  [14].  For  the  SRI,  there  may  not  be  any  beacon  light  over 
a  given  subaperture,  but  the  reference  beam  will  have  a  fraction  of  the  incoming  aperture 
averaged  intensity.  Therefore  the  local  phase  is  directly  calculated  as 


(f>  =  Tan 


-I  I' 4  -  h 


/,  -  I 


=  Tan 


1  / 1  ref  I  ref 

\  I  ref  —  I  ref 


— Tal1  (o/  0 


(2.22) 


In  this  sense,  the  SRI  is  seen  to  be  nearly  immune  to  scintillation;  however  issues  still 
arise  in  coupling  light  into  the  reference  fiber  when  these  conditions  exist.  To  minimize 
coupling  inefficiencies,  an  effective  tilt  tracking  and  correction  system  is  essential  when 
using  an  SRI.  Previous  research  also  investigated  using  doped  reference  fibers  to  achieve 
an  optical  gain  for  the  reference  beam  path.  Previous  research  by  Rhoadarmer  found  that 
just  enough  gain  to  overcome  sensor  readout  and  quantization  noise  is  ideal  [35],  but  in 
practice  this  is  not  often  implemented.  This  research  seeks  to  exploit  the  benefits  of  the  SRI 
with  respect  to  scintillation  in  conjunction  with  a  novel  wavefront  reconstruction  technique, 


discussed  later  in  Sec.  2.3. 


2.1.3  Wavefront  Reconstruction. 

It  is  often  necessary  to  further  process  the  WFS  outputs  to  achieve  an  unwrapped 
wavefront  estimate.  This  process  is  called  wavefront  reconstruction.  For  the  SRI,  WFS 
phase  outputs  are  first  sheared  into  slope  space.  Other  WFS  designs  typically  directly 
output  wavefront  gradients.  As  slopes,  each  measurement  point  is  related  to  the  others 
through  a  set  of  linear  equations  -  called  the  geometry  matrix.  This  system  of  equations  is 
then  solved  to  find  the  continuous  phase  wavefront  that  minimizes  the  squared  error  of  the 
original  system  of  equations.  The  following  subsections  describe  this  process  in  detail. 
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2. 1.3.1  Linear  Algebra. 

The  problem  of  relating  wavefront  gradient  measurements  to  unwrapped  phases  can 
be  written  as, 

Ax  =  b,  (2.23) 

where  b  is  the  measured  gradients  in  column  vector  form  of  length  m,  x  is  an  77-length 
column  vector  of  the  unwrapped  phases  being  estimated,  and  A  is  a  geometry  matrix 
relating  the  two  of  size  m  x  n.  The  specific  form  of  A  will  be  discussed  in  Sec.  2. 1.3.2, 
but  it  is  a  known  system  parameter.  The  process  of  solving  for  x  varies  depending  on  the 
specific  properties  of  A.  If  m  =  n  and  A  is  full  rank,  then  the  solution  is  straight  forward  as: 


A-1  Ax  =  A~lb 
x  =  A~xb. 


(2.24) 


However,  in  AO  this  is  rarely  the  case.  Instead,  there  are  traditionally  more  measurements 
than  phase  points  being  estimated  (m  >  /?).  This  set  of  linear  equations  is  overdetermined, 
requiring  an  approach  to  reduce  the  parameter  space  to  a  best  fit  solution.  Premultiplying  by 
A1 ,  {A1  A)  becomes  invertible,  giving  a  best  fit  solution  in  terms  of  the  minimized  squared 
error  for  x  [40]: 

ArAr  =  Arb 


(ATA)-\ATA)x  =  ( ATA)-'ATb 


x  =  (. ArA)-1Arb 

x  =  A^b 


(2.25) 


This  is  the  basic  process  which  is  used  in  the  following  subsections  to  estimate  wavefront 
phase  from  gradient  measurements  from  the  WFS.  Other  methods  exist  for  solving  a  system 
of  linear  equations  such  as  this;  however,  the  Moore-Penrose  pseudoinverse  (often  referred 
to  as  the  least  squares  pseudoinverse)  is  the  most  commonly  used  algorithm. 
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2. 1.3.2  Zonal  Reconstruction. 


Zonal  reconstruction  uses  a  simple  relative  relationship  between  measured  slopes  and 
the  phase  estimation  points.  This  relationship  is  captured  by  the  Geometry  matrix,  of  which 
there  are  two  predominant  types  depending  on  the  way  in  which  a  particular  WFS  senses 
wavefront  gradients.  The  Fried  geometry  is  typically  associated  with  the  Shack-Hartmann 
WFS,  while  the  Hudgin  geometry  is  used  with  the  lateral  shearing  interferometer  (LSI)  and 
SRI. 

Wavefront  sensors  implementing  the  Fried  geometry  sense  wavefront  slopes  at  the 
center  of  each  subaperture  and  relate  those  slopes  to  phase  estimation  points  at  the  corners 
of  each  subaperture  through  Eq.  (2.26)  and  (2.27)  [12].  Figure  2.7  shows  a  2x2  subaperture 
array  in  the  Fried  geometry. 


(.$x— Fried)  m,n  —  ^ilpm,n+ 1  (pm+ l,«+l)  ^  (fPtnji  0m+l,«) 


(^y- Fried) in ji  Pm+l,n  *Pm+\,n+  ,)“2W  m,n  4*m,n+ 1) 


(2.26) 

(2.27) 


In  this  geometry,  there  are  two  slope  measurements  per  subaperture,  creating  2 Q  equations 
for  P  phase  estimation  points.  Relating  this  to  the  discussion  in  Sec.  2. 1.3.1,  this  system  of 
linear  equations  can  be  written  as 


GF(f>  - 


t>x— Fried 


Sy— Fried 


(2.28) 


Alternatively,  wavefront  sensors  implementing  the  Hudgin  geometry  estimate  phases  at 
the  center  of  subapertures  based  on  gradient  measurements  between  adjacent  subaper¬ 
tures  [21].  Figure  2.8  demonstrates  a  2x2  subaperture  wavefront  system.  The  mathematical 
relationship  for  this  geometry  is  more  simplistic  than  the  Fried  geometry  as, 


(.$x— Hudgin)  m,n  —  4*m,n+ 1  (Pm.n  ■> 


( Sy—Hudgin)m,n  0m+ 1  ,n  * Pm,n  • 


(2.29) 

(2.30) 
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Similarly  to  the  Fried  geometry,  this  creates  2 Q  equations  for  P  phase  estimations,  creating 
the  system  of  linear  equations: 


G/[(p 


$  x- Hud  gin 


I  Sy—Hudgin\ 

Both  Eq.  (2.28)  and  (2.31)  can  be  solved  as  in  Sec.  2. 1.3.1  by: 


(2.31) 


(2.32) 


where  G  and  s  are  the  appropriate  geometry  matrix  and  slopes  from  above. 

Further,  G t  is  usually  precomputed  for  the  given  wavefront  sensing  geometry  being 
implemented.  This  eliminates  the  need  to  compute  the  pseudoinverse  at  each  frame.  The 
resulting  matrix  is  called  the  reconstruction  matrix, 

R  =  G+.  (2.33) 


This  poses  the  general  least  squares  wavefront  reconstruction  problem: 

=  Rs.  (2.34) 


2.1.4  Wavefront  Correction. 

The  process  of  correcting  the  aberrated  field  inputs  is  heavily  driven  by  the  application 
of  interest.  Techniques  exist  for  correction  of  both  field  amplitude  and  phase  [2,7,36]; 
however,  only  phase  compensation  is  discussed  further  here.  Traditional  AO  uses  the 
methods  described  in  the  previous  sections  to  estimate  the  phase  wavefronts  and  eliminate 
any  estimated  aberrations.  Figure  2.9  depicts  the  physical  process  of  phase  conjugation. 
With  an  accurate  estimation  of  the  aberrated  input,  a  conjugate  shape  is  applied  to  the 
surface  of  a  wavefront  corrector  such  that  when  the  light  reflects,  the  wavefront  is  near 
flat.  As  mentioned  previously,  the  low  order  aberrations  of  tilt  are  often  corrected 
separately  from  the  higher  order  aberrations  using  a  combination  of  deformable  mirror 
and  FSM  [25,41], 
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Figure  2.7:  The  Fried  Geometry.  Slopes  are  measured  at  the  center  of  the  subaperture  and 
are  related  to  the  phases  at  each  comer.  Slope  vectors  shown  here  include  both  x  and  y 
components  for  each  subaperture. 


There  are  a  number  of  kinds  of  wavefront  correction  devices  that  lend  to  specific 
applications  [41].  This  research  uses  a  two  inch  diameter  Optics  In  Motion  FSM  [3] 
and  12  x  12  actuator  Boston  Micromachines  Micro-Electro-Mechanical  System  (MEMS) 
DM  [1], 

2.1.5  AO  Control. 

To  implement  a  practical  AO  system,  some  form  of  control  scheme  is  necessary  to 
derive  useful  commands  to  the  wavefront  correction  device(s)  based  on  the  estimated  input 
aberrations  [25,41].  In  practice,  most  AO  systems  rely  on  closed- loop  control  for  the 
increased  robustness  when  integrating  complex  beam  control  systems.  Figure  2.10  shows 
a  high-level  control  diagram  for  a  tracking  system.  An  aberrated  beam  entering  the  AO 
system  will  have  some  time-varying  wavefront  tilt,  <pt.  The  beam  first  reflects  off  of  the 
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Figure  2.8:  The  Hudgin  Geometry.  Slopes  are  measured  between  adjacent  subapertures, 
which  are  related  to  the  phases  at  the  center  of  each  subaperture. 


FSM,  which  leaves  the  beam  with  some  residual  tilt  error,  <pt-res  (near  zero  for  a  good 
control  system  in  steady-state  operation).  The  residual  tilt  error  is  sensed  and  estimated  by 
the  tracker.  The  estimated  error  is  then  fed  into  a  control  law  to  calculate  the  next  set  of 
adjustments  to  the  FSM  -  with  the  ultimate  goal  of  driving  the  sensed  error  to  zero.  For  the 
AO  system  used  in  this  research,  a  first-order  linear  filter  of  the  form, 

ck  =  aek-i  +J3ck-U  (2.35) 

is  implemented  in  both  the  tracking  and  higher-order  control  systems.  In  Eq.  (2.35),  cy  is 
the  controller  output  at  time  step  k,  ek  is  the  sensed  error  at  time  step  k,  a  adjusts  the  DC- 
gain,  and  /?  is  a  servo  leakage  factor.  In  a  pure  integrator  setup,  (3  =  1,  but  slightly  lowering 
this  value  decreases  sensitivity  to  any  unsensed  modes  or  misregistration  errors. 

Following  the  control  law,  a  digital-to-analog  converter  (DAC)  takes  the  digitized 
controller  calculations  and  converts  them  to  a  usable  analog  voltage.  This  is  typically 
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Figure  2.9:  Wavefront  phase  conjugation  using  a  deformable  mirror. 


buffered  and/or  amplified  prior  to  being  sent  to  the  wavefront  corrector.  The  FSM  being 
used  in  this  research  implements  its  own  internal  closed-loop  controller  to  maintain  the 
mirror  position  to  the  analog  position  commands  [3].  Based  on  the  position  of  the  FSM 
with  respect  to  the  optical  axis,  a  tilt  is  imparted  to  the  incoming  beam,  and  the  process 
repeats  for  each  frame  in  time. 

2.1.6  Metrics  of  Performance. 

The  holy  grail  of  almost  any  optical  system  design  is  to  achieve  diffraction  limited 
imaging  performance  for  the  application  at  hand.  For  this  to  be  the  case,  there  must  be  zero 
wavefront  aberrations  throughout  the  optical  system  -  including  the  imaging  path.  With 
this  in  mind,  the  goal  of  AO  is  to  remove  any  observable  wavefront  aberrations  from  the 
light  of  interest  before  reaching  the  science  sensors  or  target.  In  AO,  the  primary  measure 
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Figure  2.10:  High-level  control  scheme  for  a  tilt  compensation  system. 


of  correction  performance  is  Strehl  ratio.  Strehl  ratio  is  defined  as  [38] 

0  /,>«,*/, 


(2.36) 


where  (H(fx,fy)  is  the  aberrated  optical  transfer  function  (OTF)  [18]  and  Hcu(fx,fy )  is  the 
diffraction  limited  OTF.  Both  of  these  OTF’s  are  in  terms  of  their  spatial  frequencies  fx  and 
fy.  There  are  several  approximations  to  this  form  that  better  suit  various  applications  or 
measurement  techniques  [25,38].  For  this  research  field-estimated  Strehl  is  used. 

As  the  name  suggests,  field  estimated  Strehl  takes  a  given  field  and  computes  an 
approximate  Strehl  without  the  need  to  propagate  to  a  focal  plane.  This  is  often  useful 
for  efficient  numerical  simulations  or  direct  Strehl  estimation  from  WFS  reconstructions. 
It  is  given  by 

A /]>(*, y)|2  dxdy  ) 

where  U(x,y )  is  the  field  within  an  aperture  of  area  A.  For  the  simulations  to  follow, 

Eq.  (2.37)  is  discretized  to  the  form 


\mean  [f/(x,y)]|2 
mean  [\U(x,y)\2]' 


(2.38) 


For  hardware  performance  evaluation,  field  estimated  Strehl  is  subject  to  the  accuracy  of  the 
wavefront  reconstruction.  Therefore,  when  evaluating  the  effect  of  different  reconstructors 
on  AO  performance  in  hardware,  field  estimated  Strehl  is  not  a  good  metric.  An  alternative 
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is  to  use  either  the  tracking  camera  measurements  or  a  separate,  higher  resolution,  science 
camera  to  compute  ' H(fx ,  fy)  from  the  sensed  irradiance  point  spread  functions. 

2.2  Scintillation  and  Adaptive  Optics 

In  electromagnetics,  when  a  field  interferes  with  itself,  areas  of  both  constructive  and 
destructive  interference  result.  This  is  considered  a  scintillated  field  [48].  Figures  2.1 1  and  2.12 
show  the  irradiance  of  a  point  source  beacon  propagated  through  both  weak  turbulence  and 
strong  turbulence,  respectively.  In  each  case,  D/r0  «  2  and  A  =  1.55/rm,  but  there  are  key 
differences  between  the  two  that  should  be  highlighted: 


Weak  Turbulence  Example 

•  C\  =  5  x  1(T15,  L  =  1km 

•  cr2x  =  0.01 

•  Smooth  phase  function 

•  Relatively  evenly  spread  irradiance 


Strong  Turbulence  Example 

Cl  =  2.2  x  1CT17,  L  =  244km 
cr2x  =  1.06 

Phase  function  contains  sharp  cuts 

Scintillated  irradiance  with  both  peaks 
and  areas  of  near  zero  power. 


When  describing  a  given  optical  path  through  the  atmosphere,  weak  atmospheric 
turbulence  generally  implies  that  any  transmitted  light  will  only  accumulate  a  randomized 
phase  as  a  result  of  any  refractive  or  diffractive  effects.  Figure  2.13  shows  light  being 
spherically  emitted  from  a  star  that  is  roughly  planar  upon  reaching  Earth’s  atmosphere. 
However,  once  these  planar  wavefronts  propagate  through  the  turbulent  atmosphere,  the 
wavefronts  become  aberrated.  Strong  atmospheric  turbulence  generally  refers  to  an  optical 
path  that  causes  not  only  phase  accumulation  but  also  interference  effects.  For  this  to 
happen,  the  beam  accumulates  a  phase  aberration  early  in  the  optical  path.  With  further 
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Irradiance  Unwrapped  Phase 


(a)  (b) 


Figure  2.11:  Recieved  point  source  irradiance(a)  and  phase(b)  propagated  through  weak 
turbulence. 


Irradiance  Unwrapped  Phase  (xphase) 


Figure  2.12:  Recieved  point  source  irradiance(a)  and  phase(b)  propagated  through  strong 
turbulence. 
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Figure  2.13:  Effects  of  the  randomized  atmosphere  on  planar  phase  wavefronts. 


propagation  the  aberrated  wavefronts  begin  to  constructively  and  destructively  interfere  - 
causing  scintillation  in  the  received  light. 

Unfortunately  for  AO,  scintillation  presents  a  unique  set  of  problems  with  respect  to 
phase  conjugation.  Fried  wrote  the  first  significant  publication  about  scintillation’s  effect  on 
the  corresponding  phase  function  [13],  and  he  found  that  in  areas  of  no  light  within  the  field 
a  new  phase  component  forms.  Such  points  are  called  branch  points.  Figure  2.14  shows 
an  idealized  field  that  has  a  single  branch  point  at  the  center.  At  the  branch  point  there  is 
zero  amplitude,  and  the  phase  goes  undefined.  To  make  the  field  continuous  in  a  modulo- 
one  wave  sense,  a  rotational  phase  emanates  outward  from  the  branch  point,  producing 
the  corkscrew  shape.  Expanding  to  the  general  case,  Fried  showed  that  all  branch  points 
create  a  rotational  phase  field  which  is  in  addition  to  the  generally  associated  irrotational 
wavefront  aberrations  from  propagations  through  the  turbulent  atmosphere  [13].  This  can 
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be  written  simply  as 


<f>  =  (f>irr  +  (pmt,  (2.39) 

where  <f>rot  is  the  rotational  component  of  the  phase  due  to  branch  point  effects  and  <f)lir 
is  the  irrotational  component  typically  associated  with  atmospheric  optical  propagations. 
Note  that  in  weak  turbulence,  <f>rot  =  0  because  there  are  generally  no  branch  point  effects 
present. 

As  seen  in  the  theoretical  branch  point  of  Fig.  2.14,  the  rotational  phase  is  defined 
over  a  one  wave  region.  This,  combined  with  the  rotational  nature,  requires  a  sharp  cut 
emanating  outward  from  the  branch  point  along  the  wrapping  edge  of  the  one  wave  window. 
This  is  called  a  branch  cut  and  can  be  seen  in  parts  of  the  reconstructed  phase  in  Fig.  2.12. 
Branch  cuts  always  link  pairs  of  branch  points  or  travel  from  one  branch  point  to  outside 
the  pupil  region.  The  rotational  phases  are  defined  on  an  arbitrary  one  wave  region,  and  the 
2n  window  can  be  shifted  without  changing  the  modulo-2^-  congruence  of  the  wavefront. 
However,  doing  so  moves  any  branch  cuts.  For  example,  the  rotational  phase  of  Fig.  2.14 
is  windowed  between  —n  to  n.  If  the  window  is  changed  to  0  to  In.  the  wrapping  boundary 
would  move  to  the  opposite  side  of  the  corkscrew,  and  the  branch  cut  would  shift  sides 
of  the  field.  Note  that  shifting  this  window  by  integer  multiples  of  2 n  give  equivalent 
placements  of  the  branch  cuts.  Likewise,  if  the  wrapping  window  is  kept  fixed,  a  uniform 
piston  added  to  the  entire  field  has  the  same  effect  [31,45].  This  latter  technique  is  preferred 
for  this  research  simply  due  to  computational  efficiency  and  can  be  written  as 

=  TV  [4>rot  +  h\ ,  (2.40) 

where  (f>'rot  is  still  modulo-2/r  equivalent  to  <f>rot.  Figure  2.15  gives  a  graphical  example  of 
this  process. 

Unfortunately  for  most  AO,  the  wavefront  corrector  used  is  a  continuous  face-sheet 
design.  These  correctors  will  inevitably  have  a  fitting  error  along  any  branch  cuts  in  the 
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Figure  2.14:  Theoretical  rotational  phase  function  surrounding  a  centered  branch  point  in 
a  propagated  field. 


field  being  conjugated.  The  following  section  poses  a  strategy  to  minimize  the  impact  of 
branch  cut  fitting  errors  from  the  wavefront  corrector. 

2.3  Wavefront  Reconstruction  in  Strong  Turbulence 

As  described  in  Sec.  2.1.3,  raw  WFS  measurements  usually  need  reconstruction  prior 
to  correction.  Herein  lies  a  major  problem  for  traditional  AO  techniques  when  branch 
points  are  present  in  the  observed  field.  Any  reconstruction  matrix  calculated  using  the 
least  squares  pseudoinverse  from  Eq.  (2.25)  cannot  reconstruct  a  rotational  component  [13]. 
This  is  not  a  function  of  the  particular  WFS/DM  geometry  chosen  -  it  is  instead  a  behavior 
of  the  Moore-Penrose  pseudoinverse.  Figure  2.16  demonstrates  this  behavior. 

Fortunately,  this  is  not  a  dead  end.  There  are  actually  a  number  of  algorithms  one  can 
use  to  recover  <f>rnt.  Pellizzari  provides  the  most  complete  comparison  to  date  with  a  focus 
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Figure  2.15:  Effects  of  adding  a  uniform  piston,  h,  to  the  rotational  field  while  maintaining 
a  fixed  wrapping  window.  Modulo-2^-  congruence  is  maintained,  but  the  branch  cut  moves 
as  h  changes. 


on  the  AO  application  [31].  Of  the  many  algorithms  from  a  range  of  science  applications, 
the  primary  constraint  that  must  be  considered  for  use  in  AO  is  computational  efficiency. 
The  wavefront  must  be  reconstructed  in  real-time  at  the  frame  rate  of  the  WFS  -  often  in 
the  kilo-frames  per  second.  Pellizzari  used  this  to  narrow  down  his  evaluation  to  fifteen 
reconstructor  options  for  which  he  did  a  more  thorough  comparison.  The  research  that 
follows  builds  from  his  conclusions.  In  the  subsections  below,  the  most  promising  real¬ 
time  branch  point-tolerant  wavefront  reconstructors  are  described. 

2.3.1  Post-Processing  Congruence  Operation  Reconstruction. 

The  most  promising  class  of  algorithms  is  based  on  a  post-processing  congruence 
operation  (PCO)  to  find  and  add  back  in  the  rotational  portion  of  the  field  [31,45].  This 
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Least  Squares  Reconstructed  Phase  [rad] 


Figure  2.16:  Applying  a  least  squares  reconstructor  to  a  rotational  field. 


research  chose  the  SRI  because  it  directly  senses  the  wavefront  phase  -  including  any 
rotational  components  present.  We  can  thus  formulate  the  problem  in  terms  of  Eq.  (2.39) 
and  the  PCO  technique  as 

4>irr  =  <PlS  -  RlsT  (<f>SRl)  > 

and  (2.41) 

tyrot  —  4  f 7  |  0  S  A7  in'  ]  > 

where  <pSRi  is  the  SRI’s  raw  phase  measurement  and  the  T(0)  operator  finds  the  x  and 
y  gradients  of  0,  consistent  with  the  reconstructor  geometry  being  used  [13,  15,31,45]. 
Substituting  Eq.  (2.41)  into  Eq.  (2.39)  gives 

<Ppco  -  <Pls  +  4V  [</>sri  -  (Pls]  ■  (2.42) 

This  process  maintains  the  traditional  least  squares  phase  reconstruction  technique  for  the 
irrotational  phase  component  and  finds  the  modulo-2^-  equivalent  rotational  phase  out  of 
the  original  measurement.  This  process  adds  only  a  very  small  computational  overhead  to 
gain  the  rotational  component  reconstruction. 

There  is  further  room  for  improvement  on  the  basic  PCO  process.  As  discussed  in 
Sec.  2.2,  the  branch  cut  locations  can  be  varied  by  applying  a  piston  to  the  rotational  field 
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prior  to  wrapping  Eq.  (2.40).  The  subsections  that  follow  take  the  general  PCO  algorithm 
and  attempt  to  optimize  the  particular  realization  of  rotational  phase  with  respect  to  the 
branch  cut  locations  by  seeking  an  optimal  piston  addition,  h. 

To  compare  realizations  of  rotational  phase,  a  cost  function  posed  by  Venema  [45] 
and  refined  by  Pellizzari  [31]  called  intensity  weighted  cut  length  (IWCL)  is  used. 
IWCL  first  locates  any  branch  cuts  in  a  rotational  field-again  by  shearing  the  raw  phase 
to  get  wavefront  gradients,  and  searching  for  sharp  gradients.  It  then  integrates  the 
SRI’s  estimated  field  irradiance  from  Eq.  (2.21)  along  the  estimated  cuts.  Results  from 
both  Venema  and  Pellizzari  show  high  negative  correlation  between  IWCL  and  Strehl 
ratio  [31,45].  This  research  deviates  slightly  from  the  previous  IWCL  implementations 
in  not  normalizing  the  integrated  IWCL  by  the  total  field  irradiance.  While  the  normalized 
IWCL  is  useful  for  comparing  specific  test  cases,  there  is  no  advantage  to  the  additional 
computations  in  a  real-time  AO  system. 

By  seeking  to  minimize  IWCL,  branch  cuts  are  moved  towards  areas  of  minimized 
irradiance,  thus  keeping  the  effects  of  DM  fitting  errors  along  the  cuts  to  a  minimum. 
This  is  only  one  cost  metric  -  others  can  be  chosen  based  on  the  application  of  interest. 
One  alternate  was  developed  for  beam  projection  applications,  which  is  presented  later  in 
Sec.  4.6. 

2.3.1. 1  LSPV. 


Least  squares  principle  value  (LSPV)  was  first  published  with  respect  to  AO  by 
Venema  and  Schmidt  [46].  LSPV +77  takes  the  one  wave  range  of  possible  h  offsets,  and 
divides  it  n  times: 


heval  ~ 


0  :  —  :  2/r  | . 

77 


(2.43) 


While  the  definition  above  splits  hevai  equally  across  the  one  wave  range,  there  are  variants 
that  do  not  necessarily  follow  this  constraint.  The  LSPV+/7  reconstruction  process  then 


34 


evaluates  all  n  values  of  hevai  in  terms  of  IWCL  to  determine  the  best  <f>mt: 


<f>LSPV+n  ~  &LS  +  [<f>SRI  -  <pLS  +  h<>pt\  ,  (2.44) 

hop,  =  argmin/w  [IWCL  ((f>wt(hevai))\  •  (2.45) 


LSPV+1  and  LSPV+4  were  the  primary  algorithms  investigated  in  this  level  of  PCO 
optimization.  LSPV+1  may  seem  like  a  poor  choice,  but  with  only  one  value  of  hevai,  there 
is  no  need  to  compute  IWCL  -  giving  a  computational  advantage.  LSPV+200,  or  other 
large  n,  can  be  used  to  test  for  the  true  ideal  hopt  for  a  given  field  but  is  not  suggested  for 
real-time  implementation  due  to  the  computational  requirements  and  diminishing  marginal 
gains.  Figure  2.17  shows  one  particular  LSPV+4  execution.  It  is  obvious  that  some  values 
of  hevai  are  less  optimal  than  others. 


2.3. 1.2  Iterative  Probing. 

Proposed  by  Pellizzari  [31],  this  algorithm  starts  a  narrow  LSPV+5  search  about  the 
midpoint  between  the  previous  frame’s  hopt^ i  and  ( hopt ).  In  this  case,  the  narrow  LSPV+5 
search  evaluates 


hevai  1 


-0.3  -0.15  0  0.15  0.3 


+ 


hopt.k- 1  "t"  (hopt) 
2 


waves. 


(2.46) 


The  calculated  hopUk  from  the  above  LSPV+5  search  serves  as  a  seed  for  a  second  search. 
The  new  search  narrows  the  range  of  hevai  for  a  finer  investigation  of  the  IWCL  parameter 
space  about  hoptk. 


h 


evall 


-0.07  -0.04  0  0.04  0.07 


+  hgptjc  waves. 


(2.47) 


Pellizzari  motivated  this  by  showing  that  the  IWCL-space  is  roughly  parabolic.  This 
iterative  probing  method  allows  for  further  precision  with  still  minimal  additional 
computational  requirements. 
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Figure  2.17:  These  plots  depict  the  rotational  portion  of  a  wavefront  in  radians  as  a  piston 
shift  is  applied  [45].  Branch  points  are  marked  by  X’s  and  O’s  based  on  their  respective 
charge,  and  branch  cuts  are  highlighted  with  white.  The  non-unique  nature  of  the  rotational 
phase  is  highlighted  by  applying  piston:  n  (a),  n/2  (b),  0  (c),  and  -n/2  (d).  Plots  adapted 
from  [45]. 


2.3. 1.3  Histogram  Search. 

Also  proposed  by  Pellizzari  [31],  this  algorithm  minimizes  the  area  of  <f>rot  near  the 
wrapping  boundaries.  It  first  calculates  a  histogram  of  the  <f>rot  phase  values.  From  there,  it 
calculates  h'opt  as  the  shift  requried  to  place  the  minimally  sized  histogram  bin  at  the  nearest 
wrapping  edge.  Additionally,  a  narrow  LSPV+3  search  explores  the  IWCL  parameter- 
space  about  h'opt  to  arrive  at  the  final  hopt. 
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2.3. 1.4  Brent’s  Method. 


The  final  PCO  algorithm  posed  by  Pellizzari  utilizes  what  is  known  as  a  Brent’s 
Method  [33].  Unlike  the  iterative  probing  algorithms  above,  this  type  of  search  algorithm 
minimizes  the  required  evaluations  of  the  cost  function  for  the  optimization  problem  at 
hand.  For  this,  it  means  fewer  evaluations  of  IWCL  through  a  golden  ratio  search  of  the 
h  parameter-space.  Knowing  that  the  IWCL  parameter-space  is  roughly  parabolic,  this 
algorithm  recursively  refines  a  parabolic  fit  by  adjusting  h.  Each  parabolic  fit  uses  four 
data  points,  from  which  the  search  algorithm  divides  that  search  range  in  half.  Successive 
searches  maintain  a  ratio  between  cost  function  evaluation  inputs  such  that  previous  IWCL 
evaluations  are  reused  as  the  search  narrows. 
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III.  High-Fidelity  Simulation  Development 


This  chapter  highlights  the  development  of  a  Matlab  model  of  AFIT’s  own  AO  testbed. 
From  the  start  of  development,  this  was  intended  to  be  a  high-realism  wave-optics  based 
simulation  of  the  true  system  dynamics.  The  following  sections  highlight  the  development 
decisions,  lessons  learned,  and  configuration  capabilities  for  the  simulation  model. 

The  simulations  were  developed  in  the  Matlab  programming  environment.  As  this  is 
a  discrete  numerical  model,  careful  consideration  was  given  to  adequate  spatial  sampling 
of  the  electromagnetic  fields  involved.  Although  this  discussion  does  not  cover  the  theory 
behind  numerical  simulation  of  wave-optics,  this  information  can  be  found  in  Schmidt  [38]. 
Additionally,  Matlab  toolboxes,  WaveProp  and  AOTools,  were  used  to  augment  the 
methods  proposed  by  Schmidt.  These  toolboxes  were  developed  by  the  Optical  Sciences 
Company  (tOSC).  Much  self-authored  functionality  and  integration  code  was  developed 
for  this  research  and  is  described  in  the  sections  that  follow. 

3.1  AFIT  Adaptive  Optics  System 

AFIT’s  AO  system  testbed  was  built  as  a  multi-purpose  platform  on  which  to  give  a 
hands-on  experience  for  Air  Force  students  being  introduced  to  AO.  As  it  has  evolved,  it  has 
become  an  extremely  capable  and  configurable  platform  for  advanced  AO  research  [37]. 
Despite  this,  having  a  simulation  model  on  which  to  conduct  initial  research,  where  many 
variables  of  experimental  work  are  idealized,  creates  a  powerful  time  and  money  saving 
tool.  The  following  subsections  introduce  the  physical  system  and  how  it  was  translated 
into  the  computer  model. 

3.1.1  Hardware  Specifications  and  Setup. 

Built  by  SAIC,  AFIT’s  AO  system  consists  of  two  independent  control  systems.  As 
described  in  Sec.  2.1.2,  a  tracking  system  handles  the  low  order  tilt  wavefront  aberrations, 
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while  a  separate  high-order  system  senses  and  corrects  the  remaining  wavefront  errors. 
Figure  3.1  shows  the  optical  layout  of  the  AO  system  which  can  be  fed  by  anything  that 
relays  to  an  entrance  pupil  at  the  FSM.  As  is  currently  setup,  the  beam  diameter  at  the  FSM 
is  14mm  before  being  reduced  by  a  factor  of  three  in  a  4-/  relay  to  the  DM.  From  then 
on,  the  beam  remains  at  4.67mm  and  is  relayed  once  more  to  the  sensor  entrance  pupil 
which  precedes  the  Shack-Hartman  WFS  beam-splitter.  This  size  gives  ten  DM  actuators 
across  the  beam  as  shown  in  Fig.  3.2.  A  border  of  actuators  remain  that  are  slaved  to 
the  actively  driven  ‘master’  actuators.  Slaving  unsensed  actuators  minimizes  influence  of 
adjacent  actuators  at  the  beam  edges  on  the  actively  driven  actuators  [43].  For  this  research, 
the  Shack-Hartmann  WFS  is  being  completely  bypassed  and  going  unused.  Instead,  the 
higher-order  WFS  is  a  19x19  subaperture  spatial  SRI  operating  at  1 .9Kfps  with  one  frame 
of  latency  from  sensor  exposure  to  DM  commands  out.  The  tracking  sensor  has  a  field 
of  view  of  1.8mrad  over  a  32  x  32  sensed  pixel  array  operating  at  4.76Kfps,  again  with 
roughly  one  frame  of  delay  from  exposure  to  FSM  command  out.  Note  that  this  field  of 
view  is  referenced  to  the  AO  system  entrance  pupil,  not  the  demagnified  beam  entering  the 
actual  tracker. 

This  spatially-oversampled  WFS-DM  geometry  is  important  in  properly  reconstruct¬ 
ing  the  wavefront  with  a  phase-only  WFS  when  there  is  potential  for  large  phase  differences 
between  WFS  subapertures  [6].  Such  large  phase  differences  are  a  natural  consequence 
of  branch  cuts  within  the  field.  Therefore,  a  high-resolution  wavefront  reconstruction  is 
formed  from  the  19  x  19  SRI  measurements,  which  is  then  downsampled  to  the  DM  res¬ 
olution  prior  to  the  control-law  calculations.  It  is  important  when  downsampling  only  to 
include  phase  estimates  within  the  active  aperture,  otherwise  the  outer  edge  of  downsam¬ 
pled  points  will  always  be  artificially  biased  towards  zero. 

Simulating  the  hardware  system  allows  for  a  few  simplifications.  First  of  all,  there  is 
no  need  to  simulate  any  optical  relays  of  the  entrance  pupil.  Instead,  all  phase  corrections 
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Figure  3.1:  Layout  of  AFIT’s  AO  system  testbed. 


Figure  3.2:  Geometry  of  SRI  subapertures  (dotted  squares)  to  master  DM  actuators  (blue 
dots)  the  beam  extent  is  shown  in  red. 


40 


are  performed  on  the  same  field,  which  is  then  appropriately  scaled  in  irradiance  before 
being  sensed  by  the  tracker  and  SRI.  Of  course,  this  assumes  a  perfectly  aligned  and  zero 
aberration  optical  system;  however,  if  one  wanted  to  add  either  effect,  it  would  be  a  simple 
addition. 

In  simulation,  a  one  meter  telescope  collects  light  from  the  optical  propagations.  This 
is  not  as  critical  as  D/r0  is,  where  D  is  the  entrance  aperture  diameter.  As  discussed  before, 
the  ratio  of  wavefront  corrector  size  to  Fried’s  parameter,  tq,  is  an  important  consideration 
when  designing  a  system.  For  this  simulation  work,  any  atmospheric  conditions  can  be 
generated,  leaving  the  choice  of  entrance  aperture  size  somewhat  arbitrary.  The  entrance 
telescope  is  simulated  to  have  a  focal  length  at  the  simulated  target  distance  for  each 
scenario,  giving  a  collimated  output  beam  to  the  AO  system. 

The  telescope  also  demagnifies  the  one  meter  entrance  pupil  to  a  4.67mm  exit  pupil. 
This  design  decision  ignores  the  one-third  demagnification  from  the  FSM  to  DM  in  the 
physical  system  but  increases  the  simulation  simplicity.  Due  to  the  changed  angular 
magnification  at  the  simulated  FSM,  the  only  change  between  the  simulations  and  physical 
system  is  the  final  gain  coefficient  between  the  tracker  controller  to  the  analog  FSM 
voltages.  Since  no  inertial  effects  of  the  FSM  are  simulated,  this  gain  discrepancy  has  no 
effect  on  the  fidelity  of  results.  To  justify  excluding  simulating  inertial  effects  of  the  FSM, 
a  comparison  between  the  FSM  command  frequency  spectra  (at  the  highest  Greenwood 
frequency  tested)  and  the  physical  FSM’s  measured  frequency  response  showed  only  a 
small  content  approaching  the  FSM’s  roll-off  frequency. 

With  the  optical  system  design  handled,  consideration  of  the  separate  control  systems 
must  be  paid.  Often  in  computer  simulations  of  AO  systems,  the  tracking  system  and  high- 
order  system  operate  at  integer  frequency  multiples  of  one  another  (i.e.  Tracking  at  4Kfps, 
High-order  at  2Kfps);  however,  this  is  not  often  the  case  for  physical  systems.  This  is  an 
issue  when  calculating  the  optical  fields  that  feed  the  system.  As  these  are  computationally 
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very  expensive  to  generate,  a  minimum  number  of  frames  should  be  sought.  The  solution 
used  in  this  research  was  to  calculate  separate  lists  of  sampling  times  for  each  control 
system  based  on  their  respective  sampling  rates.  With  those,  a  unique  merge  of  the  two 
lists  is  performed,  while  maintaining  a  table  that  links  each  point  in  time  to  which  sensor(s) 
is  active.  This  resulting  list  is  not  uniformly  spaced  in  time,  but  it  does  not  have  to  be. 
Knowing  the  turbulence  wind  speeds  and  point  in  time,  WaveProp  is  able  to  generate  the 
propagated  fields  without  issues.  The  listing  below  gives  pseudo  code  through  this  process 
for  additional  clarification. 
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1  %  Load  propagated  fields 

2  %  Setup  sensors  and  correctors 

3 

4  for  currFrame  =  1 :numFrames 

5  %  Apply  simulated  telescope 

6  %  Apply  last  FSM  and  DM  corrections 

7 

8  if  SRI  active 

9  %  Detect,  reconstruct,  downsample,  control  law,  filter 

10  end 

11  if  TRK  active 

12  %  Detect,  threshold,  centroid,  control  law 

13  end 

14 

15  %  Calculate  Strehl 

16  end 

With  the  above  system  coded  up,  initial  simulations  began.  Unfortunately,  the  first  set 
of  Strehl  plots  looked  like  the  left  subplot  of  Fig.  3.3.  As  the  figure  shows,  initially  the  AO 
system  is  completely  off,  and  the  Strehl  is  roughly  zero.  The  tracking  system  turns  on  0.01s, 
and  the  Strehl  picks  up  to  roughly  0.1.  At  0.02s,  the  high-order  system  turns  on,  leading 
to  a  good  initial  improvement,  but  the  system  quickly  goes  unstable.  Upon  investigation  of 
the  DM  commands,  it  was  obvious  that  tilt  was  dominating  the  DM  corrections. 

While  not  apparently  an  issue  in  the  previous  work  by  Pellizzari  where  the  separate 
control  systems  operated  at  the  same  frequencies,  here  the  higher-order  and  tilt  control 
systems  create  a  condition  where  they  fight  each  other.  As  there  is  always  going  to  be  a 
slight  residual  correction  error  from  the  tracking  system,  the  SRI  detects  this  and  tries  to 
compensate,  although  at  a  slower  sampling  rate.  The  faster  tracking  system  will  have  time 
to  sense  and  correct  for  the  same  residual  that  the  higher-order  system  is  still  processing. 
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When  the  control  loop  gains  are  high  enough,  this  behavior  can  quickly  overwhelm  the 
limited  stroke  of  the  DM  and  even  drive  the  whole  AO  system  unstable. 

The  solution  is  simple:  ensure  there  are  no  common  correction  modes  between  any 
control  loops.  In  this  case,  project  out  tilt  from  the  higher-order  system.  In  traditional  AO,  it 
is  common  to  augment  the  traditional  least- squares  reconstructor  matrix  with  a  tilt-removal 
operation,  such  that  the  reconstructed  wavefront  has  no  tilt  components.  The  process  for 
this  is  shown  in  the  code  below: 


1  %%  Remove  Z— Tilt: 

2 

3  %  Generate  arbitrary  orthogonal  tilted  wavefronts: 

4  [xt,yt]  =  meshgrid ( 1 : nDMacts ) ; 

5  %  Zero— mean: 

6  xt  =  (xt  -  mean (mean (xt ))) ; 

7  yt  =  (yt  -  mean (mean (yt ))) ; 

8  %  Only  consider  the  active  aperture: 

9  xt  =  xt (activeMask  ==  1); 

10  yt  =  yt  (activeMask  ==  1); 

11  %  Transform  from  2D— space  into  Reconstructor— space 

12  %  and  normalize  the  tilted  wavefronts: 

13  XZ  =  xt*xt '  /  (xt '  *xt)  ; 

14  YZ  =  yt*yt '  /  (yt '  *yt)  ; 

15  %  Create  the  Tilt— Removal  Operator  Matrix: 

16  P  =  eye  ( sum  (activeMask  (:))  )  -  XZ  -  YZ; 

17 

is  %  Combine  to  form  a  tilt— removed  Reconstructor: 

19  R_tr  =  P  *  R; 
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Unfortunately,  with  the  PCO-based  reconstruction  techniques,  this  process  is  not  as 
simple.  As  outlined  in  the  previous  chapter,  the  PCO  process  is 

4>PCO  =  &LS  +  [(f>SRI  -  0ls]  •  (3.1) 

However,  to  make  fipco  tilt-removed,  two  forms  of  least- squares  reconstruction  are  needed: 

<f>PCO  =  <f>LS-tr  +  [(f>SRI  -  <f>LS-ti ]  ,  (3-2) 

where  <f>Ls-tr  is  a  tilt-removed  least-squares  reconstruction  while  <f>Ls-ti  retains  tilt.  This 
does  indeed  give  a  tilt-removed  reconstruction,  but  it  requires  the  calculation  of  two  least- 
squares  reconstructions  rather  than  only  one. 

A  better  solution  exists.  Realizing  that  the  control  law  follows  the  reconstruction, 
a  design  decision  exists  on  where  to  calculate  actuator  slaving.  Often,  this  is  rolled  in 
with  the  least-squares  reconstructor;  however,  this  technique  will  not  work  with  the  PCO 
calculation  of  (f)rot.  Additionally,  the  slaving  operation  should  ideally  be  the  last  calculation 
prior  to  sending  the  DM  commands  out.  So  a  tradeoff  can  be  made  in  this  case:  only 
calculate  a  single,  tilt-included,  least-squares  reconstruction  which  gives  a  tilt-included 
PCO  reconstructed  wavefront.  Downsample  the  reconstruction,  send  it  through  the  control 
law,  and  apply  one  final  matrix  operator  which  removes  piston,  removes  tilt,  and  applies 
slaving.  This  is  the  optimal  data-flow  process  for  a  PCO  reconstructor;  shown  in  Fig.  3.4.  In 
total,  this  requires  just  two  matrix  multiplies  -  one  for  the  least-squares  reconstruction  and 
one  for  the  final  filtering.  With  the  tracking  and  higher-order  control  systems  decoupled, 
performance  improved  greatly.  This  is  seen  in  the  right  plot  of  Fig.  3.3. 

3.1.2  System  Radiometry. 

With  the  simulation  model  framework  mostly  in  place,  focus  could  be  turned  to 
modeling  radiometry  throughout  the  optical  system.  This  is  crucial  if  noise  effects  are 
going  to  be  accurately  simulated,  as  photon  shot  noise  is  signal-dependent  [17]  and  other 
noise  sources  are  only  meaningful  in  relation  to  the  signal  level. 
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Field  Estimation  Strehl 


oj  =0.04,  fg  =40Hz,  Coupled  Tilt  =0.04,  fg  =40Hz,  Decoupled  Tilt 


Figure  3.3:  The  effects  of  neglecting  tilt  decoupling  on  Strehl  ratio.  Right  plot  shows  tilt 
properly  decoupled  over  the  same  input  wavefronts  for  comparison. 


SRI 

(wrapped,  w/  tilt) 


Reconstruct  & 
Downsample 
(w/  tilt) 


Control  Law 
(low  pass  filter 
per  pixel) 


Remove 

Piston 

Remove 

Tilt 

Slave 

(Actively  illuminated  area  only) 


DAC  to  DM 


Figure  3.4:  Optimal  process  flow  for  a  PCO  based  reconstruction  algorithm. 
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3. 1.2.1  Radiometric  Assumptions. 

In  the  physical  AO  system,  power  can  be  varied  to  each  signal  path  using  polarized 
beam  splitters  preceded  by  adjustable  half-wave-plates  [37].  This  requires  the  input  light 
be  initially  filtered  by  a  polarization  analyzer.  For  the  simulation  work,  no  polarization 
effects  are  modeled  -  only  the  power  splitting  ratios.  In  this  sense,  the  input  light  for  the 
simulations  would  be  equivalent  to  just  after  the  initial  polarization  analyzer  on  the  physical 
table. 

A  number  of  assumptions  were  made  in  this  radiometric  model: 

1.  No  losses  due  to  reflections. 

2.  No  losses  due  to  dispersion  or  misalignments. 

3.  Ideal  beam  splitters. 

Realistically  there  is  a  small  loss  through  the  system  despite  the  anti-reflective  coatings 
used  on  all  optics.  Additionally,  the  beam  splitters  have  a  finite  extinction  ratio.  This  is 
relevant  for  the  beam  paths  that  theoretically  receive  no  light.  Realistically,  the  Shack- 
Hartmann  and  auxiliary  beam  paths  receive  a  small  percentage  of  the  light. 

Additionally,  realistic  sensor  models  were  created  in  WaveProp  based  on  the  camera 
datasheets  from  Goodrich.  For  both  the  tracker  and  SRI,  a  Goodrich  SU320-KTSW  sensor 
is  used,  with  the  key  specifications  shown  in  Table  3.1. 

3. 1.2.2  Realistic  SRI  Reference  Fiber  Coupling. 

Previous  simulations  assumed  the  SRI  had  a  constant  fiber  coupling  coefficient  into 
its  single-mode  reference  fiber,  regardless  of  the  input  wavefront  condition.  This  is  a 
rough  approximation  for  quick  evaluation  but  is  insufficient  for  a  realistic  noise  analysis, 
determining  low  signal  performance,  and  examining  transient  behavior  before  the  AO 
system  reaches  steady-state. 
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Table  3.1:  Goodrich  SU320-KTSW  Sensor  Specifications 


Resolution 

320  x  256t 

Pixel  Pitch 

25/rm 

Quantum  Efficiency 

>0.65tt 

Full  Well  e~ 

750,000e~ 

RMS  Read  Noise 

183e-ttt 

I  Can  address  a  smaller  rectangular  region  of  the  array  for  increased  speed. 

II  Actual  quantum  efficiency  at  1550nm  is  likely  greater,  but  0.65  was  used  for  these 
simulations.  A  sensor  calibration  can  experimentally  determine  this  value  [22]. 

"  '  Based  on  correspondence  with  Goodrich  technicians  for  a  FPA  gain  of  one. 

Based  on  previous  work  by  Wheeler  [47],  frame-by-frame  fiber  coupling  efficiency 
and  exit  mode  shape  can  be  calculated.  Using  the  physical  optic  and  fiber  parameters  of 
AFIT’s  SRI  hardware,  the  single  mode  accepted  into  the  fiber  (LP0i)  is  back-propagated  to 
the  SRI  entrance  pupil.  With  this,  the  input  aberrated  beam  is  projected  onto  the  known 
back-propagated  LP0i  fiber  mode,  giving  a  realistic  frame-by-frame  coupling  coefficient. 
Similarly,  the  exiting  LP0i  fiber  mode  is  a  magnified  version  of  the  input  mode,  based  on 
the  specific  re-collimating  optics  used  after  the  output  of  the  fiber.  This  output  mode  is 
scaled  by  the  current  frame’s  coupling  coefficient  as 

Ufibemutik )  —  (Uinput(k),  Uq \  -in  )U0 1  -out  9  (3.3) 

where  (•,  •)  denotes  an  inner  product,  k  represents  the  current  time  step,  Ujnput  is  the  input 
aberrated  beam  of  the  SRI  reference  leg  at  time  step  k,  f/0 i_,-„  is  the  normalized  back- 
propagated  LPqi  mode  that  is  coupled  into  the  fiber,  and  f/0i_OMf  is  a  normalized  version 
of  the  collimated  fiber  exit  mode.  This  operation  captures  the  stochastic  variations  in 
coupling  efficiency,  while  other  effects  like  Fresnel  reflections  are  treated  separately  and 
deterministically.  As  is  shown  in  Figs.  3.5  and  3.6,  fiber  coupling  is  particularly  important 
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in  modeling  transient  behavior  immediately  after  closing  the  loop  in  strong  turbulent,  low 
signal  situations. 

3. 1.2. 3  Power  Splitting. 

As  shot  noise  is  dependent  on  the  signal  level  at  the  detector  and  the  other  noise 
sources  are  only  meaningful  relative  to  the  signal  level,  careful  work  had  to  be  taken  to 
ensure  proper  radiometry  throughout  the  simulated  AO  system.  This  problem  was  initially 
approached  by  adjusting  the  power  splitting  between  tracker  and  SRI  to  achieve  roughly 
equivalent  signal  levels  on  each  detector.  For  the  tracker  this  is  straight  forward,  but  the 
SRI  presents  a  challenge  with  the  input-dependent  fiber  coupling.  For  this  calibration, 
a  low-Rytov  turbulence  scenario  was  used  to  determine  the  optimum  tradeoff  in  power 
splitting.  Figures  3.5  and  3.6  illustrate  this  tradeoff.  Before  the  AO  system  is  engaged, 
the  fiber  coupling  is  near  zero.  First  the  tilt  system  is  engaged,  giving  a  very  slight 
increase  in  coupled  light.  This  coupled  signal  must  be  sufficient  to  form  the  SRI  wavefront 
estimate,  otherwise  the  higher-order  control  loop  may  not  successfully  close  when  engaged. 
When  the  higher-order  system  engages,  coupling  immediately  jumps  as  the  DM  starts 
correcting.  This  is  shown  in  Fig.  3.5  when  the  fiber  power  becomes  stronger  than  the 
signal  power.  However,  as  fiber  coupling  efficiency  is  directly  correlated  with  Strehl  ratio, 
stronger  turbulence  degrades  the  power  coupled  into  the  SRI  reference  fiber,  as  seen  in 
the  medium-Rytov  case  of  Fig.  3.6.  A  balance  must  be  struck  between  supplying  enough 
power  to  initially  engage  the  higher-order  control  loop  while  not  taking  too  much  power 
from  the  signal  leg.  The  power  splitting  can  thus  be  optimized  for  the  expected  atmospheric 
conditions,  or  as  in  this  research,  a  general  compromise  was  found  that  accommodates  a 
wide  range  of  scenarios.  Figure  3.7  shows  the  power  splitting  breakdown  through  the 
AO  system  that  was  empirically  determined,  and  as  it  turns  out  matches  up  closely  to  the 
physical  AFIT  hardware  setup. 
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Signal  vs  Fiber  Output  Power  Splitting  [refSplit  =  0.8]  Fiber  Coupling  Efficiency 


Figure  3.5:  SRI  reference  fiber  coupling  as  the  AO  system  starts  up.  The  left  plot  shows 
the  power  in  the  fields  for  both  the  signal  and  reference  (labeled  Fiber  Power)  legs  through 
the  SRI.  The  right  plot  calculates  the  fiber  coupling  efficiency  over  time.  Atmospheric 
conditions  are  fg  =  20Hz,  Rytov  number  =  0.04. 


Signal  vs  Fiber  Output  Power  Splitting  [refSplit  =  0.8]  Fiber  Coupling  Efficiency 


Figure  3.6:  SRI  power  splitting  between  reference  and  signal  legs  on  the  left  using  a 
LSPV+4  reconstructor  for  turbulence  with  a  Rytov  number  of  0.488  and  Greenwood 
frequency  of  20Hz.  The  resulting  fiber  coupling  is  shown  on  the  right. 
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Figure  3.7:  Approximate  power  splitting  within  the  simulated  AO  system.  The  SRI’s 
interferograms  are  represented  with  I\  through  /4.  I,  and  Is  are  the  reference  and  signal 
paths,  respectively.  The  fiber  loss  will  depend  on  the  instantaneous  Strehl  ratio  of  the  field 
entering  the  SRI. 
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3.1.3  Sensor  Noise  Effects. 

With  the  system  radiometry  setup,  focus  can  be  turned  to  the  noise  sources.  For  this 
discussion,  focus  is  on  a  single  pixel  of  one  SRI  interferogram,  as  we  are  most  interested 
in  noise  effects  relating  to  the  reconstruction  performance.  However,  it  should  be  noted 
that  the  balanced  power  distributions  and  same  detector  models  used  in  the  SRI  and  tracker 
yield  roughly  equivalent  SNR’s.  Noise  sources  analyzed  include  photon  shot  noise  cr]hot, 
detector  readout  noise  cr2read,  and  ADC  quantization  noise  cr2quant^ation  so  that  the  total  system 
noise  is  given  by 


0~sys  — 


+  cr 


read 


+  CT 


2 

quantization  * 


(3.4) 


Photon  shot  noise  follows  a  Poisson  distribution  and  is  only  a  function  of  the  signal 
strength.  Its  noise  contribution  to  Eq.  (3.4)  can  be  expressed  as  [35] 


Odiot  =  K 

rice  [  ,r  m  (Is  +  TJcIr)  ]  ^ 

~  M  1  in,J^PX  r. 

robins  L 

where  K  is  the  mean  photoelectron  count  for  a  given  pixel,  rjqe  is  the  detector  quantum 
efficiency,  Nhms  is  the  number  of  SRI  interferograms  captured  simultaneously,  Tint  is  the 
detector  integration  time  in  seconds,  3Kpx  is  the  area  of  one  pixel  in  square  meters,  rjc  is  the 
mean  coupling  efficiency  into  the  single-mode  fiber,  and  Is  and  lr  are  the  irradiances  for  the 
signal  and  reference  legs  of  the  SRI,  respectively.  As  this  is  an  interferogram,  Is  and  /,.  are 
averaged,  rather  than  simply  adding  their  signals. 

As  mentioned  in  Table  3.1, 


thread  =  083e~)2  • 


(3.6) 


This  read  noise  number  varies  depending  on  the  gain  and  calibrations  used,  but  this  gives 
a  rough  middle  ground  for  an  InGaAs  camera  in  this  class  and  represents  one  bit  of  digital 
resolution.  This  sensor  uses  a  12bit  ADC  to  quantize  the  analog  photocurrents.  Assuming 
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a  rounding  ADC,  the  quantization  noise  can  be  expressed  as  [35] 


cr 


2 

quantization 


T2LSB2 

^(18302 


(3.7) 


where  there  are  183  photoelectrons  per  least  significant  bit  (LSB)  of  the  12bit  ADC. 

With  all  of  the  noise  sources  defined,  the  challenge  becomes  to  define  a  meaningful 
SNR  to  test  AO  and  reconstructor  performance.  The  noise  is  now  defined  in  terms  of 
photoelectrons  for  an  individual  pixel  of  one  interferogram.  The  signal  will  thus  be  defined 
as  the  average  number  of  photoelectrons  incident  on  one  pixel  over  the  same  integration 
period.  Making  use  of  Fig.  3.7,  roughly  5%  of  the  entrance  power  reaches  the  pixel  from 
both  the  signal  and  reference  legs  of  the  SRI.  Since  this  is  an  interference  pattern,  the 
combined  signal  could  range  between  0-10%  with  a  mean  about  5%.  Therefore,  the  signal 
estimate  is  5%  of  the  input  power.  Note  that  this  decreases  if  Strehl  performance  degrades 
for  any  reason.  Putting  Eqs.  (3.4)-(3.7)  together  with  the  previous  definition  of  signal,  SNR 
can  be  defined  as 


SNR 


SNR  = 


K 


0-2  _|_  0-2 

shot  read  quantization 

K 


jK  +  (me-)2  +  (^me-f 


(3.8) 

(3.9) 


Using  the  above  equation  for  SNR,  the  required  input  beam  power  can  be  calculated 
for  1 ,55/jm  light.  Note  that  these  SNR  values  are  not  in  decibels.  In  the  extreme,  before 
the  higher-order  loop  is  closed,  actual  signal  levels  can  be  much  lower  than  these  estimates 
due  to  poor  fiber  coupling,  as  seen  in  Figs.  3.5  and  3.6.  Because  of  this,  initially  closing  the 
loop  has  been  found  to  be  the  primary  obstacle  for  the  SRI-based  AO  system  to  overcome. 
Ensuring  the  tracking  system  is  closed  before  trying  to  close  the  higher-order  loop  is 
essential  in  low  signal  conditions. 
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3.2  Atmospheric  Propagations 


While  not  the  primary  focus  of  this  research,  realistic  aberrated  optical  inputs  are 
needed  to  test  the  AO  system  over.  An  ideal  point  source  beacon  at  a  wavelength  of  1 .55//m 
was  used  in  all  cases.  As  mentioned  before,  care  was  taken  in  ensuring  adequate  spatial 
sampling  of  the  wavefronts  using  methods  from  Schmidt  [38].  WaveProp  also  alerts  users  if 
inadequate  sampling  is  detected.  Table  3.2  lists  the  spatial  sampling  used  in  this  research, 
where  cr~  is  the  Rytov  number,  d/r0  is  the  ratio  of  the  DM  actuator  spacing  to  Fried’s 
parameter,  6src  is  the  beacon  source  plane  grid  spacing,  SAosys  is  the  receive  plane  grid 
spacing  entering  the  simulated  telescope,  N  Grid  Points  is  the  number  of  discrete  spatial 
samples  per  side  of  the  simulated  field,  C\  is  the  index  of  refraction  structure  constant  along 
the  propagation  path,  and  D-  is  the  propagation  path  length. 

3.2.1  Parameter  Exploration. 

To  fully  characterize  the  AO  system  performance,  a  wide  parameter-space  was 
investigated.  Most  importantly,  this  research  evaluates  the  effect  of  an  increasing  Rytov 
number  on  reconstruction  performance.  Additionally,  for  each  Rytov  number,  a  range  of 
Greenwood  frequencies  and  SNR’s  are  investigated.  Finally,  two  r0’s  were  considered  at 
d/ro  =  1  (Challenging)  and  J/r0  =  1/2  (More  realistic).  This  is  all  summarized  in  Fig.  3.8. 


Table  3.2:  Atmospheric  propagation  setup  parameters. 


d/r0 

6src  [mm] 

dAosys  [mm] 

N  Grid  Points 

c2 

D-  [km] 

0.04 

1 

0.4 

8 

1024 

2e-14 

1 

0.5 

1 

4 

4 

2048 

le-15 

20 

1 

1 

5 

7 

2048 

4e-16 

50 

0.5 

0.5 

9 

13 

1024 

7.88e-17 

80 

1 

0.5 

12 

15 

2048 

2.2e-17 

244 
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Horizontal  Propagation  of  a  point  source 
Hold  d/r„  constant  in  all  cases:  1,  V2 
Vary  Rytov  Number:  0.04,  0.5,  1 
Vary  Greenwood  Frequency:  20,  40,  60,  80,  100Hz 
Lead  to  various  combinations  of  Dz,  C%,  and  Vl 


wind 


rj*-f 


-fMr 

W 


r  2 


C 2 


Kvind 


D, 


Figure  3.8:  Point  source  atmospheric  propagation  setup. 
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IV.  Simulation  Results 


This  chapter  first  presents  a  compiled  set  of  simulation  results  from  the  various 
reconstruction  methods  introduced  in  Ch.  II.  Obvious  performance  improvements  from  the 
traditional  least  squares  reconstruction  approach  are  seen,  especially  as  turbulence  strength 
increases.  While  that  was  the  original  intended  scope  for  the  simulation  portion  of  this 
research,  further  investigations  led  the  author  to  several  additional  insights  on  PCO-based 
reconstruction  that  greatly  clarify  previous  questions  and  hypotheses.  These  new  findings 
include  least  squares  phase  anomalies  due  to  gradient  aliasing,  the  cause  and  potential 
solutions  to  sharp  fades  when  using  PCO  reconstruction  in  strong  turbulence,  and  refined 
behavior  models  of  hopt  during  closed-loop  operation.  This  chapter  concludes  with  an 
extension  of  the  PCO  technique  in  strong  turbulence  to  beam  projection  applications. 

4.1  Reconstructor  Performance 

Due  to  the  large  number  of  simulations  run,  this  section  presents  a  summarized  set 
of  the  raw  results.  While  exploring  the  effects  of  SNR  and  Greenwood  frequency  were 
informative  along  the  research  pathway,  the  more  general  performance  trends  between 
each  reconstruction  technique  can  be  explored  without  these  details.  The  tables  that 
follow  present  averaged  statistics  of  the  field-estimated  Strehl  ratio  from  each  one  second 
simulation  run.  For  each  Rytov  number  and  reconstructor  combination,  all  simulations 
were  grouped  together  to  calculate  the  mean,  standard  deviation,  and  normalized  standard 
deviation  of  the  field-estimated  Strehl  ratio.  The  mean  Strehl  ratio  measures  absolute 
AO  correction  performance,  and  the  standard  deviation  of  the  Strehl  ratio  measures  the 
stability  of  the  AO  correction.  Higher  mean  Strehl  ratios  and  lower  standard  deviations  are 
desired.  More  detailed  plots  of  the  simulations  run,  linking  effects  of  SNR  and  Greenwood 
frequency,  are  found  in  Appendix  A. 
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Table  4.1:  Averaged  field-estimated  Strehl  ratios  over  all  Greenwood  frequencies  and 
SNR’s  tested  for  the  cases  listed. 


LS 

LSPV+1 

LSPV+4 

Histogram 

Iter.  Probe 

Brent’s 

0.04 

0.5406 

0.5416 

0.5416 

0.5416 

0.5416 

0.5416 

0.5 

0.3242 

0.3844 

0.3869 

0.3876 

0.3872 

0.3846 

1 

0.0949 

0.2047 

0.2177 

0.2183 

0.2202 

0.2134 

Table  4.2:  Field-estimated  Strehl  ratio  standard  deviation,  averaged  over  all  Greenwood 
frequencies  and  SNR’s  tested  for  the  cases  listed. 


LS 

LSPV+1 

LSPV+4 

Histogram 

Iter.  Probe 

Brent’s 

0.04 

0.0329 

0.0322 

0.0322 

0.0322 

0.0322 

0.0322 

0.5 

0.0835 

0.0498 

0.0496 

0.0470 

0.0487 

0.0488 

1 

0.0677 

0.0647 

0.0674 

0.0625 

0.0626 

0.0646 

Table  4.3:  Field-estimated  Strehl  ratio  normalized  standard  deviation,  averaged  over  all 
Greenwood  frequencies  and  SNR’s  tested  for  the  cases  listed. 


LS 

LSPV+1 

LSPV+4 

Histogram 

Iter.  Probe 

Brent’s 

0.04 

0.0628 

0.0613 

0.0613 

0.0613 

0.0613 

0.0613 

0.5 

0.2581 

0.1310 

0.1295 

0.1223 

0.1267 

0.1278 

1 

0.7162 

0.3198 

0.3152 

0.2899 

0.2881 

0.3058 
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4.2  Low  SNR  Reconstruction  and  Phase  Discrepancies 

Among  the  first  observations  made  was  in  the  initial  weak  turbulence  testing.  In 
weak  turbulence,  there  are  no  branch  point  effects  present,  and  therefore  (f)rot  =  0  in  the 
absence  of  any  noise  effects.  The  expectation  was  to  observe  equivalent  performance 
between  all  reconstructors.  However,  a  consistent  performance  improvement  of  all  tested 
PCO  reconstructors  was  seen  over  the  least  squares  reconstructor,  specifically  at  the  lowest 
SNR’s  tested.  This  was  very  much  a  surprise,  as  prior  work  hypothesized  that  the  addition 
of  noise  may  be  debilitating  to  the  PCO  technique  [31].  This  sparked  a  series  of  follow-up 
investigations  to  explore  what  was  going  on.  All  tests  performed  verified  the  initial  findings 
that  all  the  PCO  reconstructors  outperformed  the  least  squares  reconstructors  in  low-SNR 
situations.  The  question  then  became:  Why? 

An  initial  theory  was  proposed  by  the  author  that  was  based  on  differing  noise 
propagation  through  the  two  reconstruction  techniques.  The  hypothesis  is  based  on  the 
idea  that  the  least  squares  operator  could  actually  amplify  the  SRI  measurement  noise  such 
that 

\\<t>Lsn\\2  >  m\2  (4.D 

where  (f>LSlI  is  the  least  squares  propagated  noise  component,  (f>„  is  the  noise  component 
from  the  measurement  process,  and  ||  •  H2  is  the  Frobenius  norm  [16].  This  idea  was  formed 
by  calculating  the  condition  number  [16]  of  the  least  squares  reconstruction  matrix  and 
finding  it  to  be  roughly  15  -  ignoring  the  suppressed  reconstruction  modes.  Now  relate  this 
to  the  PCO  algorithm  defined  in  Eq.  (2.42): 

<f>pco  =  &LS  +  '4V  [</>sri  -  <!>ls]  •  (4.2) 

The  above  form  can  be  expanded  in  terms  of  all  noise  components  involved: 

<t>pco  -  fiLSj,  +  0ls„  +  '4V  {TV  [(f)  +  (f>„]  -  ((f>LSt  +  0z.s„)}  •  (4-3) 
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The  theory  proposes  that  the  two  <f>LSn  terms  cancel  each  other  out  in  the  PCO  process, 
leaving  only  (f>„  contributing  to  the  final  reconstruction.  Under  the  assumption  of  Eq.  (4. 1), 
only  having  (f>„  would  be  less  of  a  degradation  than  A  rough  analytical  trace  for 

this  theory  was  created,  but  it  needed  additional  verification  to  justify  some  simplifications 
made. 

A  simulation  was  written  to  model  this  process  which  proved  the  theory  wrong.  This 
verification  simulation  models  AO  system  behavior  in  steady-state  closed-loop  operation 
by  generating  an  irrotational  test  field  (could  be  as  simple  as  a  plane  wave)  that  has 
a  nearly  flat  wavefront.  The  mean  field  irradiance  is  scaled  using  the  equations  from 
Sec.  3.1.3  to  have  an  SNR  of  roughly  1.5.  This  field  is  sent  through  the  SRI  model  to 
add  realistic  measurement  noise  effects,  and  that  field,  <f>SRi,  is  sent  through  both  a  least 
squares  reconstructor  and  a  simple  LSPV+1  reconstructor.  From  these,  the  noise  terms  can 
be  extracted  and  examined. 

Specifically,  the  simulation  found  that 

</>ls„  -  <f>n,  (4.4) 

which  invalidated  the  author’s  original  theory  for  the  PCO  reconstructor  performance 
improvement.  In  hopes  of  gleaning  more,  additional  realizations  of  noise  were  run, 
and  something  unexpected  was  noted  in  a  number  of  cases.  Amongst  other  things,  the 
simulation  plotted  <f>LSn  -  (/>,,.  usually  showing  no  difference.  However,  Figs.  4.1a  and  4.1b 
show  two  examples  of  what  appeared  to  be  branch  points  at  first  glance. 

Yet,  these  are  not  branch  points  for  a  number  of  reasons.  First,  the  field  that  was  input 
to  the  SRI  had  a  uniform  irradiance.  One  could  argue  that  shot  noise  creates  areas  of  zero 
irradiance,  which  is  a  valid  argument.  The  detected  field  irradiance  can  be  calculated  from 
the  SRI  interferograms  by  Eq.  (2.21),  and  Fig.  4.2  shows  the  estimated  input  irradiance 
corresponding  to  the  measurement  of  Fig.  4.1b.  It  is  clear  that  areas  of  zero  detected 
photons  are  scattered  all  over  the  aperture,  but  there  is  no  correlation  to  the  anomalies 
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observed  in  Fig.  4.1b.  Since  these  effects  do  not  seem  to  be  caused  by  the  input  irradiance, 
looking  to  the  SRI  phase  measurement,  a  clue  is  found.  Shown  in  Figs.  4.3a  and  4.3b,  the 
points  from  Figs.  4.1a  and  4.1b  can  be  traced  to  points  in  equivalent  locations  within  the 
SRI  phase  that  make  transitions  larger  than  n  radians  between  adjacent  subapertures  due 
to  measurement  noise.  While  examining  <(>sri,  closer  examination  was  given  to  both  the 
least  squares  and  LSPV+1  reconstructed  wavefronts,  and  it  was  clear  that  while  the  PCO 
reconstructor  successfully  reconstructed  the  original  phase,  the  least  squares  reconstruction 
had  anomalies  at  the  locations  of  the  sharp  phase  differences  in  <Psri- 

To  further  examine  this  phenomenon,  a  second  simulation  was  run  that  isolates  the 
noise  generated  large  phase  differences  from  the  SRI.  A  hypothetical  (pSRi  is  generated  that 
is  completely  flat  except  for  two  adjacent  center  subapertures  that  are  set  to  -0.9n  and  0.9;r, 
as  shown  in  Fig.  4.4a.  Without  the  surrounding  additional  noise,  this  simulation  exposes  the 
aperture-wide  reconstruction  effects  due  to  large  noise-induced  adjacent  phase  difference. 
Figure  4.4c  shows  the  LSPV+1  reconstruction,  matching  the  SRI  wavefront  with  an  extra 
n  piston  shift  from  the  PCO  h  term.  On  the  other  hand,  the  least  squares  reconstruction 
clearly  identifies  the  source  of  the  phase  anomalies  evident  in  Fig.  4.1. 

The  (f>rot  term  from  the  PCO  algorithm  cancels  out  the  least  squares  anomaly;  however, 
this  exposes  a  misnomer  in  this  case.  The  wrapped  term  in  the  PCO  algorithm,  which  in  a 
noiseless  case  is  equivalent  to  any  rotational  component  in  the  phase  function,  can  actually 
be  any  phase  components  missed,  or  in  this  case  artificially  introduced  by  the  least  squares 
reconstruction  relative  the  SRI  phase  measurements.  The  author  likens  this  to  accounting 
for  slope  discrepancies  when  using  a  Shack-ffartmann  WFS.  Since  the  SRI  measures  phase, 
not  wavefront  gradients,  the  author  terms  this  as  accounting  for  the  phase  discrepancy.  In 
the  case  of  the  noise  induced  phase  discrepancies,  they  are  not  actually  rotational  in  nature. 
As  described  in  Sec.  2.3  a  least  squares  generated  operator  cannot  reconstruct  rotational 
phases.  Further,  Fig.  4.5  plots  the  gradients  for  one  of  these  anomalies,  where  red  represents 
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a  positive  field  and  blue  is  negative.  The  author  stresses  these  arguments  to  say  these  are 
not  so-called  “noise-induced  branch  points”.  They  are  not  rotational,  and  they  are  actually 
anomalies  introduced  by  the  least  squares  reconstruction  that  do  not  exist  in  the  true  phase 
function. 

Since  these  phase  discrepancies  result  from  noise  creating  artificially  large  adjacent 
phase  measurements  on  an  otherwise  nearly-flat  wavefront,  it  is  now  clear  why  all  PCO 
algorithms  outperform  the  least  squares  reconstructor  when  no  branch  points  are  present 
and  only  at  the  very  lowest  SNR’s.  The  author’s  original  theory  failed  to  explain  why  the 
performance  gain  was  only  seen  at  SNR’s  lower  than  five,  as  it  predicted  a  lowered  noise 
propagations  for  any  noise  present.  The  new  theory  requires  such  low  SNR’s  to  induce 
the  large  adjacent  phase  measurement  errors.  If  only  using  a  least  squares  reconstructor  in 
these  cases,  the  control  system  will  attempt  to  conjugate  these  phase  discrepancies  that  do 
not  exist  in  the  true  wavefront  leading  to  potentially  unstable  AO  performance. 

The  cause  for  the  observed  least  squares  behavior  turns  out  to  be  an  aliasing  effect 
of  the  wavefront  gradients  fed  to  the  the  least  squares  reconstructor.  Section  4.5  presents 
the  research  that  led  to  this  recent  understanding.  Any  gradients  larger  than  n  radians  are 
assumed  to  be  wrapping  cuts;  however  these  points  are  not  wrapping  cuts.  So,  the  large 
adjacent  gradients  are  wrapped,  but  the  surrounding  gradients  remain.  This  is  the  cause  for 
the  phase  discrepancy  surrounding  such  points. 


4.3  Closed-Loop  Behavior  of  hopt 

Pellizzari’s  work  found  a  numerical  PDF  approximation  for  hopt  [31],  which  found 
(hopt)  =  0.  However,  there  is  a  very  important  difference  between  preceding  research  and 
this  work  in  the  implementation  of  the  PCO  wrapping  operator.  Previously,  the  wrapping 
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(a)  (b) 


Figure  4.1:  Different  realizations  of  unexpected  differences  between  LS  noise  term  and 
PCO  noise  term. 


SRI  Irradiance  Estimate 
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Figure  4.2:  Irradiance  estimate  of  input  SRI  field  from  noise  behavior  simulation.  This 
corresponds  to  the  measurement  of  Fig.  4.1b. 
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Figure  4.3:  SRI  phase  measurements  of  a  tilted  wavefront  where  SNR«  1.5.  In  the 
subplots,  (a)  corresponds  to  Fig.  4.1a  and  (b)  corresponds  to  Fig.  4.1b. 


was  performed  by 


W[0]  = 


Tan  1  {imag  [exp  (i2n<f>)\ ,  real  [exp  (ilncp)]] 
2n 


(4.5) 


which  is  physically  correct  but  not  computationally  efficient.  Switching  to  a  modulus 
operator,  this  process  can  be  implemented  easily  on  most  CPU’s  and  FPGA’s  used  for 
AO  control.  This  is  simply  written  as 


TV  [0]  =  mod  (0,  2 n) . 


(4.6) 


Both  methods  are  correct,  but  the  wrapping  window  for  Eq.  (4.6)  is  now  from  0-2^  radians, 
rather  than  -n-n  radians.  When  this  research  began  testing  LSPV  + 1 ,  based  on  Eq.  (4.6),  the 
AO  performance  was  terrible.  After  recalculating  the  PDF  of  hopt,  it  was  obvious  that  the 
PDF  calculated  by  Pellizzari  needs  to  be  defined  with  respect  to  the  wrapping  window  being 
used.  By  shifting  his  entire  PDF  to  be  over  the  new  wrapping  window,  (hopt)  =  n.  It  turns 
out  the  initial  implementation  of  FSPV+1  tested  here,  used  the  least  optimal  h,  leading 
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Figure  4.4:  Simulation  of  low-SNR  noise  effects  between  least  squares  and  LSPV+1 
reconstructors. 


Figure  4.5:  Simplified  wavefront  gradient  diagram  of  a  noised  induced  phase  anomaly.  The 
center  two  squares  represent  the  central  discrete  pixels  of  the  anomaly. 
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to  the  poor  results  observed.  This  is  an  important  lesson  for  practically  implementing 
these  algorithms  or  anyone  implementing  a  non-standard  wrapping  function  for  a  PCO 
application. 

Over  the  course  of  this  research,  it  was  noted  that  the  LSPV+1  reconstructor 
demonstrated  near  equivalent  AO  performance  to  the  more  advanced  PCO  reconstructors 
tested.  The  computational  simplicity  of  LSPV+1  makes  it  a  promising  lead  for  high-speed 
and  low-latency  AO.  It  was  also  noted  that  the  LSPV+1  performance  gap  between  the  other 
PCO  algorithms  widens  as  the  Rytov  number  increases.  Although  the  performance  is  still 
very  similar,  as  seen  by  Table  4.3,  the  increasing  performance  gap  between  LSPV+1  and 
the  other  PCO  reconstructors  at  a2  =  1  prompted  an  investigation  as  to  why. 

LSPV+1  always  uses  h  =  (hopt),  which  in  the  weak  Rytov  case  is  more  than 
sufficient  when  no  branch  points  are  present.  However,  as  the  Rytov  number  increases  and 
additional  rotational  phase  components  are  introduced,  one  can  surmise  that  only  having 
one  realization  of  (f)rot  could  be  limiting  in  some  cases.  The  previous  research  did  not 
show  significant  variations  in  the  PDF  of  hopt  as  the  Rytov  number  increased  above  weak 
turbulence,  perhaps  due  to  the  absence  of  noise  effects.  However,  a  series  of  simulations 
were  run  to  retest  this  idea  with  the  hypothesis  that  the  PDF  of  hopt  broadened  away  from 
(hopt)  as  the  Rytov  number  increased.  To  generate  the  PDF’s,  LSPV+200  was  run  over 
a  series  of  one-second  cases,  each  with  increasing  challenge  for  the  AO  system.  For 
each  simulation,  a  histogram  was  generated  of  the  hopt  values  chosen  by  the  LSPV+200 
algorithm  to  generate  a  rough  PDF.  Figures  4.6  and  4.7  show  the  simulation  results  with 
the  histograms  as  well  as  field  estimated  Strehl  ratio  plots  overlaid  with  the  normalized 
hopt  as  a  function  of  time.  As  the  AO  beacon  conditions  worsen,  the  histogram  shows  a 
broadening  of  hopt' s  PDF.  All  turbulence  parameters  showed  a  slight  broadening  effect  of 
the  distribution,  but  increasing  Rytov  number  is  by  far  the  dominant  factor. 
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The  confirmed  broadening  of  the  hopt  PDF  explains  why  the  LSPV+l-LSPV+4 
performance  gap  widens  as  the  atmospheric  conditions  worsen.  This  fact  also  may  lend 
to  choosing  a  specific  PCO  reconstructor  based  on  expected  conditions.  As  the  conditions 
become  increasingly  difficult,  particularly  with  respect  to  Rytov  number,  it  may  be  more 
beneficial  to  take  the  additional  computational  hit  with  one  of  the  other  PCO  reconstructors. 

4.4  Fades  in  Strong  Turbulence 

While  debugging  the  reconstruction  and  AO  system  simulation,  a  mechanism  was 
implemented  to  record  videos  of  many  important  system  parameters  as  the  system  stepped 
through  time.  Following  the  development  phase,  this  capability  was  used  to  study  the 
overall  behavior  of  the  AO  system  when  using  a  PCO  wavefront  reconstructor.  This 
additional  work  resulted  in  identifying  the  source  of  fades  when  using  a  PCO  reconstructor 
when  a  significant  rotational  phase  is  present  in  the  input  wavefronts. 

At  system  startup,  the  WFS  observes  the  fully  aberrated  beacon  wavefronts,  (p.  At 
this  point,  the  complete  rotational  phase  from  any  branch  points  reaches  the  WFS.  When 
the  AO  system  switches  on,  the  DM  conjugates  the  reconstructed  WFS  field.  Now,  the 
WFS  senses  only  the  residual  wavefront  error,  (presid  -  not  the  fully  aberrated  wavefronts. 
The  control  law  maintains  the  DM  commands  at  roughly  the  conjugate  of  the  incoming 
aberrations,  but  the  wavefront  reconstructor  is  no  longer  seeing  the  fully  aberrated  field. 
As  (presid  flattens  out,  the  PDF  of  hop,  becomes  increasingly  centered  about  (hopt). 

This  was  all  known  prior  to  the  research  presented  here;  however,  additional  insights 
may  be  gleaned  by  thinking  deeper  into  the  process.  In  steady-state,  closed-loop  operation, 
the  residual  wavefront  error  is  near-zero.  Any  branch  points  and  corresponding  branch 
cuts  that  are  present  within  (p  may  not  be  visible  in  (preside  and  even  if  they  are  partially 
apparent,  they  would  most  likely  be  of  reduced  visibility  due  to  partial  compensation.  In 
this  regime,  the  branch  cut  finding  algorithm  may  easily  miss  cuts  that  are  actually  present 
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Figure  4.6:  LSPV+200  Strehl  ratio  and  hopt  plots  over  one  second  simulation  runs  to 
evaluate  the  behavior  of  hopt' s  PDF  as  atmospheric  conditions  change. 
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Figure  4.7:  LSPV+200  Strehl  ratio  and  hopt  plots  over  one  second  simulation  runs  to 
evaluate  the  behavior  of  hopt' s  PDF  as  atmospheric  conditions  change. 
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in  the  corrected  field.  The  author  poses  a  question  as  to  what  IWCL  really  means,  since 
IWCL  was  meant  to  be  calculated  over  the  fully  aberrated  wavefront,  (f).  not  <f>resid- 

So  what  is  going  on  when  a  PCO  algorithm  changes  hopt  while  in  closed-loop 
operation?  Well,  it  depends  primarily  on  the  rotational  effects  present  in  the  field;  in  other 
words,  the  Rytov  number.  Starting  the  discussion  with  a  weak  Rytov  number  case,  the 
rotational  component  of  the  PCO  reconstruction  was  seen  to  be  almost  always  completely 
zero,  with  the  exception  of  any  noise  induced  effects  presented  in  Sec.  4.2.  In  this  weak 
turbulence  case,  hopt  very  rarely  changes,  and  when  it  does  change  due  to  noise  induced 
phase  anomalies,  it  has  no  effect  on  the  overall  system  performance. 

In  the  medium  Rytov  number  cases  simulated,  similar  behavior  is  noted.  Figure  4.8 
shows  a  plot  of  field-estimated  Strehl  ratio,  overlaid  with  each  frame’s  hopt  for  an  LSPV+4 
reconstructor.  In  this  regime,  hopt  changes  frequently,  but  overall  AO  performance  remains 
unaffected.  This  initially  diagreed  with  the  author’s  theories,  as  the  PCO  rotational 
phase  component  was  seen  to  have  consistent  notable  content.  However,  the  most  recent 
discoveries  show  that  this  content  is  not  rotational  in  nature.  The  cause  is  discussed  at 
length  in  Sec.  4.5.  Without  a  significant  rotational  component  within  the  residual  phase 
error,  hopt  varies  in  a  way  to  keep  the  PCO  phase  discrepancies  within  the  wrapping  window 
and  away  from  the  wrapping  boundaries.  At  a  change  in  hopt,  only  the  reconstructed 
wavefront’s  piston  changes,  but  since  piston  is  removed  prior  to  sending  DM  commands 
out,  there  is  no  performance  impact.  While  there  are  a  few  exceptions,  this  is  the  primary 
behavior  observed  in  the  medium  Rytov  number  cases. 

Things  change  with  strong  Rytov  number  cases.  As  Fig.  4.9  shows,  any  changes  in 
hopr  mostly  result  in  sharp  fades.  More  specifically,  the  fades  distinctly  follow  changes  in 
hopt.  Using  the  video  capability  developed  for  debugging,  it  was  noted  that  the  tracking  spot 
shape  is  not  effected  at  these  fades  but  is  shifted  off-axis.  This  behavior  was  consistently 
observed  at  the  fades  in  Strehl  ratio.  It  seemed  that  the  change  in  hopt  somehow  propagates 
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Figure  4.8:  The  top  plot  shows  field  estimated  Strehl  ratio  and  hopt  over  time  with  an 
LSPV+4  reconstructor.  The  bottom  plot  shows  the  IWCL  for  each  of  the  four  4>rot 
realizations  over  the  same  time  interval.  This  is  in  conditions  of  cr^  =  0.5  ,fg  =40Hz, 
d/r0  =  1/2. 
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Figure  4.9:  The  top  plot  shows  field  estimated  Strehl  ratio  and  hopt  over  time  with  an 
LSPV+4  reconstructor.  The  bottom  plot  shows  the  IWCL  for  each  of  the  four  4>rot 
realizations  over  the  same  time  interval.  This  is  in  conditions  of  cr2x  =  1 ,  f„  =40Hz, 
d/r0  =  1/2. 


a  change  in  tilt  through  the  AO  system,  thus  causing  the  fades.  This  prompted  several 
simulation  tests,  one  of  which  lead  to  an  explanation. 

A  number  of  AO  simulations  were  rerun  with  the  LSPV+4  reconstructor,  this  time 
calculating  the  tilt  components  for  each  PCO  rotational  phase  realization.  Since  hopt 
seemed  to  be  at  the  root  of  the  behavior,  the  least  squares  reconstruction  was  ignored  for 
these  tilt  calculations.  Figure  4.10  shows  the  results  from  the  three  Rytov  cases.  First  to 


71 


note  is  that  the  tilt  components  of  the  PCO  rotational  phase  definitely  differ  as  h  is  varied! 
For  the  weak  case,  all  rotational  tilt  components  vanish  as  soon  as  the  higher-order  control 
loop  closes  -  as  expected.  In  the  medium  case,  the  three  most  likely  /-z-valucs  produce 
small,  yet  similar  tilts.  However,  in  the  strong  case,  even  after  the  higher-order  control  loop 
engages,  there  is  widely  varying  tilt  components  between  each  rotational  phase  realization. 
This  is  the  key. 

In  the  strong  Rytov  case,  when  the  PCO  algorithm  changes  hopt,  a  new  wavefront 
tilt  realization  is  reconstructed.  When  this  propagates  through  the  rest  of  the  higher-order 
control  data  path,  described  in  Sec.  3.1.1,  the  final  tilt-removal  operator  cuts  out  whatever 
tilt  component  exists.  To  the  tracking  system,  these  changing  tilts  are  certainly  measurable. 
This  is  the  cause  for  the  fades  in  Strehl  ratio.  While  the  higher-order  correction  of  the 
wavefronts  do  not  suffer  during  the  changes  in  hopt,  it  takes  the  tracking  system  roughly  two 
frames  to  account  for  the  tilt  changes  from  the  high-order  system.  To  explicitly  demonstrate 
this  process,  Appendix  B  traces  a  hypothetical  rotational  phase  through  the  AO  system. 

4.4.1  Potential  Fixes. 

All  the  tilt  components  are  known  variables,  and  as  such,  a  solution  to  prevent  the 
fades  is  certainly  possible.  In  fact,  there  are  likely  a  number  of  potential  solutions.  One 
method  would  be  to  implement  a  feed-forward  tilt-offloading  path  to  the  tracker  from  the 
higher-order  system.  Another  method  is  to  allow  the  higher-order  control  system  to  correct 
for  some  or  all  of  the  reconstructed  tilt.  At  first,  the  author  recalled  the  tilt  coupling  issues 
presented  in  Sec.  3.1.1,  but  by  lowering  the  control  system  gains,  a  stable  state  can  be 
reached.  While  this  imposes  larger  stroke  requirements  on  the  deformable  mirror,  it  is  a 
relatively  simple  controls  problem  compared  to  the  feed-forward  implementation  with  two 
different  operating  frequencies. 

The  AFIT  AO  system  model  was  modified  to  compare  AO  performance  over 
equivalent  realizations  of  high  Rytov  turbulence  with  both  a  tilt  coupled  and  decoupled 
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cr^  =  0.04,  d/ro  =  1,  SNR=60,  fg  =  40Hz 


Time  [s] 


<7^  =1,  d/ro  =  1,  SNR=60,  =  40Hz 


Figure  4.10:  Tilt  magnitudes  for  each  PCO  rotational  phase  realization  of  a  LSPV+4 
reconstructor  over  time  for  weak,  medium,  and  strong  Rytov  numbers.  The  higher-order 
control  system  is  seen  to  engage  at  t  =0.008s. 
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control  system.  The  limitations  of  such  a  setup  are  made  clear  by  selecting  three 
increasingly  aggressive  control  gains.  Figure  4.11  shows  the  results,  where  the  left  column 
is  decoupled  and  the  right  column  is  coupled.  The  gain  starts  at  relatively  conservative 
points,  giving  good  performance  in  each  case.  However,  note  there  are  no  fades  in  the 
coupled  case,  even  as  hopt  changes!  Increasing  the  gain,  the  decoupled  performance 
sharpens  but  the  coupled  performance  begins  to  go  unstable.  Further  increasing  the  gain, 
the  decoupled  system  still  runs  well,  but  the  coupled  system  quickly  drives  itself  unstable. 
Note  that  this  test  in  simulation  allows  for  perfect  component  alignments.  In  reality  there 
will  always  be  a  slight  misalignment  between  components  and  sensors.  Because  of  this, 
a  zero  DC  gain  controller  would  be  required  on  both  the  tracker  and  higher-order  control 
system  to  implement  a  coupled  tilt  control  setup. 

This  test  verifies  the  theory  on  tilt  propagation  from  the  changing  PCO  rotational  phase 
realizations  to  the  tracking  system  as  the  cause  for  the  sharp  fades  previously  observed  with 
these  reconstructors.  Further,  suggestions  for  alleviating  this  are  given. 

4.5  High  SNR  Phase  Discrepancies 

A  very  recent  modification  to  the  simulations  was  to  separate  and  plot  (f>LS  and  <prot 
as  the  AO  frames  increment.  To  the  author’s  great  surprise,  phase  discrepancies  similar 
to  those  described  in  Sec.  4.2  were  clearly  apparent  in  high  SNR  conditions  of  all  Rytov 
numbers  tested.  This  section  describes  the  author’s  initial  tests  and  theories  as  to  why 
these  phase  discrepancies  exist  and  what  impact  they  have  on  AO  performance.  The  author 
stresses  that  these  are  only  initial  theories,  and  further  work  is  needed  to  refine  these  ideas. 

Figure  4.12  shows  an  example  of  the  AO  system  in  weak  turbulence,  with  an  SNR 
of  60,  prior  to  closing  any  of  the  correction  loops.  In  such  low  Rytov  number  and  high 
SNR  conditions  the  PCO  rotational  phase  term  was  expected  to  have  no  content,  but  that 
was  not  the  case.  The  first  clue  came  when  the  control  loops  started  correcting  in  the  weak 
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TR:  rej  =1.047,  J,  =G0Hz,  recon  1,  SNR-60,  19  SRI  Apa,  Ttk|0.99,  44]  SHI|0.997,  SO] 


rej  =1.047,  /,  =60Hz,  recon  1,  SNR-60,  19  SRI  Apa,  Ttk|0.99,  44]  SHI|0.997,  80] 


TR:  oj  =1.047,  f,  =G0Hz,  recon  1,  SNR-60,  19  SRI  Apa,  Ttk[0.99,  44]  SRI[0.097,  90] 


<rj  =1.047,  f ,  =60Hz,  recon  1,  SNR-60,  19  SHI  Apa,  Ttk[0.99,  44]  SRI[0.997,  90] 


TR:  oj  =1.047,  /,  =G0Hz,  recon  1,  SNR-60,  19  SRI  Apa,  Ttk[0.9»,  44]  SRI[0.997,  100]  =1.047,  f,  =60Hz,  recon  1,  SNR-60,  19  SRI  Apa,  TVk|0.99,  44]  SRI|0.997,  100] 


Figure  4.11:  Plots  of  AO  performance  with  tilt  decoupled  (left  column)  and  tilt  coupled 
(right  column)  between  the  two  control  systems.  Each  row  implements  an  increasing 
higher-order  control  gain. 
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turbulence  case.  Within  several  frames,  the  PCO  rotational  component  became  flat  -  as  we 
had  expected  for  the  entire  low-Rytov  simulation,  shown  in  Fig  4.13. 

The  primary  differences  with  regard  to  the  reconstruction  process  between  the  two 
above  cases  are  the  magnitudes  of  the  residual  wavefront  gradients.  These  phase  effects 
are  also  present  in  the  higher  Rytov  number  cases,  but  the  low  Rytov  case  excludes  any 
rotational  field  effects  from  the  analysis.  A  number  of  simplified  tests  were  conducted 
to  look  at  how  high  localized  wavefront  gradients  propagate  through  the  reconstruction 
process.  The  first  question  was  to  explore  how  the  SRI  and  reconstructors  handle  slopes  of 
increasing  magnitude.  As  previously  described,  the  slopes  that  the  reconstructors  operate 
on  are  generated  by  shearing  the  raw  SRI  measurements  in  the  x  and  y  directions.  To 
account  for  wrapping  cuts,  the  sheared  slopes  are  wrapped  to  be  within  -n  <  s  <  n  radians. 
This  slope  wrapping  operation  sets  up  the  fundamental  limit  on  slope  magnitudes  that 
can  be  properly  reconstructed  by  the  least  squares  process.  Any  slopes  between  adjacent 
subapertures  larger  than  n  radians  will  be  interpreted,  incorrectly,  as  a  wrapping  cut. 

Figure  4.14  explores  this  idea  with  two  tilted  wavefronts  in  the  first  row  of  plots. 
The  corresponding  SRI  measurements  of  each  tilted  field  are  shown  in  the  middle  row, 
and  the  final  row  shows  each  least  squares  reconstruction.  The  first  column  is  correctly 
reconstructed,  but  the  second  column,  which  exceeds  the  tilt  sampling  limit,  has  two 
significant  errors.  First,  the  magnitude  of  tilt  over  the  wavefront  is  greatly  reduced  from  its 
actual  tilt.  Second,  the  direction  of  tilt  is  reversed.  The  author  hypothesizes  that  the  degree 
of  each  of  these  errors  will  change  based  on  how  much  tilt  is  being  passed  into  the  system 
-  similar  to  aliasing  in  discrete  signal  theory  [29]. 

Figure  4.15  extends  the  previous  aliasing  test  closer  to  the  atmospheric  simulation 
conditions.  A  flat  wavefront  with  a  sharp,  linear  slope  down  the  center  of  the  field 
is  created.  The  slope  is  at  roughly  0.87t  radians  per-SRI  subaperture.  Additionally, 
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Figure  4.12:  Plots  of  AO  state  in  conditions  of  cr2  =  0.04,  d/r0  =  1.  Both  tracking  and 


high-order  control  loops  are  open. 


77 


Received  Beacon  Irradiancce 


Received  Beacon  Phase 

180 

120 

60 

0 

-60 

-120 

-180 

-2  -1.5  -1  -0.5  0  0.5  1  1.5  2  (deg) 


(mm) 

Residual  Phase  Error 


■* 

.  M 

l|  ^  TjjTi  YT 
„  .  #  .  ■ 

V 

■jhjL 

ry 

► 

#  +  J*  & 

dP*1 

r 

Hr 

fc* 

pi 

120 


-60 


(deg) 


LS  Reconstructed  Phase 


SRI  Phase  Measurement 


1-120 

ri  m 


i 


120 


PCO  Rotational  Phase  Term 


U 


(deg) 


-120 


.  1 

™40 

-20 

r  || 

-0 

-20 

■  -  -40 

1 

180 


120 


(deg) 


(deg) 


Figure  4.13:  Plots  of  AO  state  in  conditions  of  cr2  =  0.04,  d/r0  =  1.  Both  tracking  and 
high-order  control  loops  are  closed. 
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Figure  4.14:  Inducing  aliasing  with  the  SRI  and  least  squares  reconstructor.  The  left 
column  shows  roughly  0.97T  radians  of  tilt  between  adjacent  subapertures,  while  the  right 
column  has  roughly  1.5;r  radians  of  tilt  between  adjacent  subapertures. 
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a  central  region  of  the  slope  is  set  at  roughly  1.5/r  radians  per  subaperture.  This  is 
similar  to  most  of  the  atmospheric  wavefront  being  within  the  slope  sampling  limit  of 
the  reconstructor,  with  small  regions  that  exceed  it.  The  first  observation  made  was  that 
the  least  squares  reconstruction  mostly  resembles  the  input  wavefront,  but  it  is  clear  there 
are  some  discrepancies,  including  aliasing  over  the  high-slope  region.  In  the  initial  tests, 
only  the  LSPV+1  algorithm  was  used  for  comparison.  While  it  does  produce  a  correct 
modulo-2^-  congruent  reconstruction,  the  (f>ro,  realization  created  has  sharp  cuts  that  are  less 
than  ideal  for  any  continuous  face-sheet  DM.  This  observation  motivated  a  test  with  the 
LSPV+4  reconstructor,  which  calculated  the  most  accurate  reconstruction.  To  the  PCO 
algorithms,  the  sharp  cuts  from  the  LSPV+1  <f>rot  realization  are  treated  as  if  they  were 
branch  cuts.  In  the  LSPV+4  reconstructor,  the  cuts  are  penalized,  giving  a  realization  that 
minimizes  the  sharp  transitions. 

The  author  claims  these  tests  point  out  that  sharp,  localized  wavefront  gradients  are 
the  cause  for  the  observed  phase  discrepancies  when  reconstructing  a  high-SNR  field.  The 
LSPV+4  rotational  phase  term  from  Fig.  4.15  has  a  high  resemblance  to  the  features  seen 
in  the  atmospheric  simulations.  While  these  features  were  not  found  in  low-Rytov  number 
closed  loop  operation,  they  were  found  in  the  higher  Rytov  number  closed  loop  simulations. 
Figures  4.16  and  4.17  show  closed  loop  operation  for  cr2y  =  0.5  and  cr2  =  \,  respectively. 
The  strong  turbulence  plots  show  the  system  as  it  is  initially  closing  the  loop,  and  some 
rotational  phase  due  to  branch  points  is  likely  present  in  addition  to  the  phase  anomalies 
being  discussed  here. 

These  tests  shed  additional  light  on  the  increased  performance  drop  between  LSPV+1 
and  the  other  PCO  reconstructors  as  the  Rytov  number  reaches  one.  The  broadening 
of  hopt  s  PDF  must  be,  in  part,  due  to  these  phase  discrepancies.  Further,  the  previous 
simulations  performed  by  Pellizzari  would  likely  not  have  encountered  these  issues  with  a 
simulated  SRI  of  over  twice  the  subaperture  density.  While  increasing  the  SRI  subaperture 
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Figure  4.15:  Inducing  aliasing  with  the  SRI  and  least  squares  reconstructor.  A  tilt  runs 
down  the  center  of  the  wavefront,  and  the  very  central  portion  exceeds  the  SRI  tilt  sampling 
limit.  This  is  traced  through  the  SRI  and  various  reconstructors. 
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density  is  a  potential  solution  to  this  behavior,  one  must  consider  the  additional  splitting 
of  light  and  its  impact  on  SNR.  The  initial  results  shown  here  demonstrate  an  unintended 
benefit  of  the  PCO  algorithm  as  it  compensates  for  reconstruction  errors  due  to  localized 
wavefront  gradient  aliasing. 

4.6  Application  to  Beam  Projection 

All  research  up  to  this  point  regarding  PCO-based  wavefront  reconstructors  has  been 
studied  with  respect  to  one-way  AO  systems.  Such  one-way  applications  include,  but  are 
not  limited  to,  remote  sensing  and  astronomy.  Another  class  of  applications  involve  dual¬ 
path  propagation,  namely  those  which  not  only  receive  light,  but  also  project  compensated 
outgoing  beams.  Beam  projection  applications  include  optical  communication  and  laser 
weapon  systems.  These  applications  typically  involve  more  horizontal  propagation  paths, 
thus  entering  the  strong  turbulence  regime  much  more  frequently  than  vertical  path 
applications.  This  section  presents  initial  results  and  observations  for  a  range  of  wavefront 
correction  techniques  with  respect  to  beam  projection  applications. 

The  author  hypothesizes  most  of  the  same  conclusions  made  for  the  single-path  AO 
problem  hold  true  but  introduces  an  alternative  to  intensity  weighted  cut  length  (IWCL) 
for  PCO  reconstruction.  The  author  hypothesizes  that  cut  length  (CL)  is  a  better  PCO 
reconstruction  metric  for  beam  projection  applications  where  the  outgoing  beam  is  not 
amplitude  modulated  to  match  the  received  beacon  field.  While  all  PCO  reconstructed 
fields  will  be  modulo-2^-  equivalent,  different  PCO  variations  differ  in  the  placement  of 
detected  branch  cuts  within  the  field.  Assuming  an  ideal  corrector,  this  does  not  actually 
make  a  difference.  However,  when  a  continuous  face-sheet  DM  is  used,  there  will  always 
be  a  DM  fitting  error  in  trying  to  reproduce  the  sharp  branch  cuts. 

For  initial  comparison  of  wavefront  reconstruction  techniques  with  respect  to  beam 
projection  applications,  a  simplified  simulation  was  developed  using  tOSC’s  WaveProp. 
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Figure  4.16:  AO  system  state  during  closed  loop  operation  in  conditions  of  cr2  =  0.5, 
d/r0  =  1,  SNR=60.  LSPV+4  is  the  wavefront  reconstructor. 
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Figure  4.17:  AO  system  state  as  the  higher-order  control  loop  is  closing  in  conditions  of 


cr~x  —  1,  d/r0  =  1,  SNR=60.  LSPV+4  is  the  wavefront  reconstructor. 
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Again,  an  ideal  point  source  beacon  illuminates  the  turbulent  path  of  interest,  with  a  Rytov 
number  of  one.  The  receive  telescope  is  40  centimeters  in  diameter,  focused  at  the  beacon 
distance.  An  idealistic  SRI  senses  the  collimated  wavefront  (no  fiber  coupling  or  noise 
effects),  and  an  ideal  DM  applies  the  appropriate  conjugation  to  the  outgoing  beam.  The 
ideal  DM  has  no  stroke  limit,  negating  the  need  for  a  FSM  or  tracking  system.  No  temporal 
effects  are  included  in  this  initial  experiment.  Instead,  unique  realizations  of  turbulent  paths 
are  generated,  through  which  the  beacon  propagates.  For  each  turbulent  realization,  the  AO 
system  senses  the  field  and  various  reconstruction  techniques  calculate  corrections  for  the 
outgoing  beam.  Each  of  these  pre-compensated  beams  are  then  back-propagated  through 
the  turbulence  to  the  beacon  source  plane,  where  the  irradiance  patterns  are  averaged  over 
all  turbulence  realizations  for  each  reconstruction  method.  The  author  acknowledges  the 
lack  of  temporal  effects  will  not  paint  a  complete  picture  of  expected  AO  performance. 
However,  this  is  an  important  first  step  to  both  develop  a  dual-path  test  simulation,  as  well 
as  compare  general  reconstructor  performance  in  these  difficult  atmospheric  conditions. 

The  following  tests  compare  reconstruction  performance  of  traditional  least  squares, 
LSPV+20  using  IWCL  as  a  metric,  LSPV+20  using  CL  as  a  metric,  and  Xphase  as 
developed  by  tOSC.  LSPV+20  was  chosen  instead  of  LSPV+4  in  order  to  more  accurately 
compare  the  PCO  technique  with  respect  to  each  metric  function.  For  reference,  ideal  phase 
conjugation  and  ideal  field  conjugation  are  also  tested.  The  first  round  of  simulation  results 
is  shown  in  Fig.  4.18.  To  the  author’s  initial  surprise  there  was  no  apparent  difference 
between  the  ideal  phase  conjugation,  PCO  methods,  and  Xphase.  This  points  out  an 
important  point:  these  reconstruction  algorithms  all  produce  modulo-2^-  equivalent  results, 
but  each  uses  a  different  rotational  realization  in  hopes  of  optimizing  placement  of  any 
branch  cuts  with  respect  to  an  expected  DM  fitting  error.  This  initial  simulation  did  not 
include  any  DM  inter-actuator  coupling  effects.  Further,  the  SRI  and  DM  were  assumed  to 
have  the  same  spatial  resolution  as  the  simulated  fields  themselves.  Without  any  of  these 
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realistic  effects,  all  of  the  above  reconstruction  methods  gave  exactly  equivalent  results. 
The  poor  performance  of  the  least  squares  reconstruction  and  great  performance  of  the 
field  conjugations  also  serve  as  good  sanity  checks  before  moving  forward. 

The  simulation  was  modified  to  include  realistic  SRI  sampling  and  DM  parameters 
(except  actuator  stroke).  As  this  introduces  some  subjectivity  for  the  parameters  chosen, 
a  range  was  evaluated  to  explore  the  current  realistic  parameter-space.  Inter-actuator 
couplings  of  both  10%  and  20%  were  evaluated  for  DM’s  of  size  10  x  10,  19  x  19,  and 
32  x  32.  The  respective  SRI  subaperture  arrays  were  19  x  19,  37  x  37,  and  63  x  63, 
again,  following  the  sampling  recommendation  set  by  Barchers  [6].  The  WFS  to  DM 
command  downsampling  method  was  a  matrix  operator  based  on  a  pyramid  convolution 
three  subapertures  wide,  only  considering  the  active  subapertures.  The  ideal  phase  and 
field  conjugation  methods  remain  at  the  full  field  sampling  density  for  comparison. 

Again,  100  realizations  of  turbulence  were  run  for  each  of  the  six  DM  setups  described 
above.  The  results  of  which  are  shown  in  Figs.  4.19-4.24.  Additionally,  Figs.  4.25  and  4.26 
provide  a  visual  example  of  the  ‘best’  and  ‘worst’  DM  fits  for  an  example  turbulence 
realization  for  each  reconstructor.  These  initial  results  show  some  interesting  behaviors; 
although,  not  all  good.  Most  notably,  the  10  x  10  DM  seems  to  outperform  the  higher 
density  DM’s.  Upon  reinspection  of  the  simulation  code,  the  AFIT  DM  setup  was  used  - 
including  a  guardband  ring  of  DM  actuators  which  is  not  used  in  the  higher  density  DM 
setups.  While  it  is  not  definitive  that  the  guardband  was  the  cause  for  the  apparent  improved 
performance,  it  is  a  difference  that  must  not  be  included  in  future  tests.  Secondly,  the 
LSPV  algorithms  and  Xphase  often  show  higher  peak  on-axis  power  than  the  ideal  phase 
conjugation.  The  author  hypothesizes  this  is  a  result  of  a  constructive  interference  effect; 
although,  that  is  only  a  theory.  One  final  concern  upon  viewing  these  results  was  the  varying 
power  levels  between  each  simulation  case.  In  the  initial  code,  the  turbulence  seeds  were 
set  based  on  the  computer  clock  -  resulting  in  the  observed  behavior.  For  future  work,  the 
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seeds  are  now  set  based  on  the  trial  number  such  that  each  DM  simulation  will  use  the  same 


set  of  turbulence  realizations. 

These  initial  results  show  the  field  conjugation  performing  best  and  least  squares 
worst.  Beyond  that,  Xphase  seems  to  slightly  edge  out  LSPV+20  IWCL  in  most  cases. 
Surprisingly,  the  LSPV+20  using  CL  as  the  reconstruction  metric  underperformed  the 
IWCL  metric.  The  author  hypothesizes  this  is  another  interference  effect  at  play,  but  is 
uncertain  on  the  precise  means.  However,  before  making  too  many  conclusions  based  on 
the  results  at  hand,  the  author  stresses  the  flaws  pointed  out  in  the  preceding  paragraph. 
Further  simulations  are  being  conducted  to  form  a  better  result-set  for  comparison. 
Additionally,  future  simulations  including  temporal  effects  are  crucial  to  better  understand 
this  problem. 
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Figure  4.18:  Target-plane  irradiance  patterns,  averaged  over  100  turbulence  realizations 
with  an  ideal  DM.  The  central  circles  show  the  diffraction  limited  spot  size. 
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Figure  4.19:  Target-plane  irradiance  patterns  using  a  10  x  10  DM  with  10%  inter-actuator 
coupling. 
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Ideal  Field  Conjugation 
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Figure  4.20:  Target-plane  irradiance  patterns  using  a  19  x  19  DM  with  10%  inter-actuator 
coupling. 


90 


Ideal  Field  Conjugation 
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Figure  4.21:  Target-plane  irradiance  patterns  using  a  32  x  32  DM  with  10%  inter-actuator 
coupling. 
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Ideal  Field  Conjugation 


Ideal  Phase  Conjugation 
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Figure  4.22:  Target-plane  irradiance  patterns  using  a  10  x  10  DM  with  20%  inter-actuator 
coupling. 
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Figure  4.23:  Target-plane  irradiance  patterns  using  a  19  x  19  DM  with  20%  inter-actuator 
coupling. 
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Figure  4.24:  Target-plane  irradiance  patterns  using  a  32  x  32  DM  with  20%  inter-actuator 
coupling. 
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Ideal  Field  Conjugation 
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Figure  4.25:  Example  DM  profiles  using  a  10  x  10  DM  with  20%  inter-actuator  coupling. 
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Figure  4.26:  Example  DM  profiles  using  a  32  x  32  DM  with  10%  inter-actuator  coupling. 
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V.  Conclusions 


This  chapter  presents  a  summary  of  the  research  completed  as  part  of  this  Master’s 
work.  Challenges  overcome,  key  results,  and  future  work  are  given.  The  research  and 
observations  made  here  build  upon  the  previous  work  of  Pellizzari  and  Schmidt,  where  this 
effort  was  aimed  at  transitioning  their  research  towards  a  realistic  implementation.  In  doing 
so,  many  new  discoveries  were  made  in  the  hopes  of  extending  current  AO  capabilities  in 
strong  turbulence  regimes. 

5.1  Challenges  Overcome 

•  Developed  realistic  AO  system  model.  This  was  the  crucial  next  step  from  the 
previous  research.  A  realistic  SRI  model  was  developed,  incorporating  frame-by- 
frame  fiber  coupling  effects,  a  first  of  its  kind  for  this  research.  System  radiometry 
was  modeled,  allowing  for  meaningful  shot  noise,  senor  read  noise,  and  quantization 
noise  effects.  Care  was  given  to  correctly  model  the  independent  tilt  and  higher- 
order  control  systems,  including  the  complete  higher-order  data  processing  pipeline. 
In  doing  so,  data-flow  optimizations  were  made  specific  to  PCO-based  wavefront 
reconstruction  and  AO  control. 

•  Complex,  multi-dimensional  parameter  space.  This  research  analyzed  effects 
of  Rytov  number,  Fried’s  parameter,  Greenwood  frequency,  and  SNR  for  each 
wavefront  reconstructor  tested.  This  five  dimensional  parameter  space  created 
hundreds  of  gigabytes  of  simulation  results  which  had  to  be  efficiently  tracked  and 
managed.  Running  the  simulations  to  explore  the  complete  parameter  space  took 
weeks,  even  utilizing  several  high  performance  computers.  Given  the  simulation 
requirements,  significant  development  effort  was  given  to  optimize  code.  This 
framework  is  now  available  for  future  research. 
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•  Bring  AF1T  AO  hardware  to  an  operational  state.  As  this  research  began,  AFIT’s 
AO  system  was  plagued  with  a  variety  of  issues  that  led  to  unreliable  performance. 
Working  in  conjunction  with  the  primary  contractor,  the  author  assisted  in  finding, 
diagnosing,  and  correcting  a  number  of  software  and  hardware  issues.  The  author 
learned  Lab  VIEW  in  order  to  become  proficient  with  the  hardware  interfacing  and 
control  code.  Since  no  manuals  were  available,  the  author  learned  and  wrote  quick 
guides  for  common  procedures.  A  full-system  optical  alignment  procedure  was 
developed,  reducing  the  time  required  from  12  hours  to  only  two.  Unfortunately, 
at  the  final  stages  of  this  research  a  circuit  board  failure  rendered  the  DM  inoperable 
before  any  meaningful  data  could  be  taken.  Initial  AO  system  characterization  results 
were  presented  at  the  2011  SPIE  Defense,  Sensing,  and  Security  conference. 

•  Developed  an  ATS  configuration  change  to  allow  high-Rytov  simulations.  In  the 
configuration  from  the  contractor,  the  ATS  could  simulate  high-Rytov  number 
conditions  but  only  with  prohibitively  small  r0’s.  In  the  simulations  performed  the 
ratio  of  DM  actuator  spacing  to  r0  was  kept,  at  most,  to  one.  However,  in  the  original 
configuration,  the  ATS  could  not  maintain  that  ratio  at  cr~  =  1.  The  author  developed 
an  optical  modification  that  would  allow  much  greater  emulation  capabilities  in  the 
r0  :  cr2  parameter  space  of  interest  to  this  research. 

5.2  Key  Results 

•  Demonstrated  greatly  improved  strong  turbulence  AO  performance  in  terms  of  both 
mean  Strehl  and  Strehl  variance.  The  AO  performance  gains  using  a  PCO  wavefront 
reconstructor  over  the  traditional  least  squares  reconstructor  in  terms  of  mean  field 
estimated  Strehl  ratio  and  Strehl  ratio  normalized  standard  deviation  are  summarized 
in  Table  5.1.  Mean  Strehl  is  a  metric  of  absolute  AO  correction  performance  (higher 
is  better),  and  the  Strehl  standard  deviation  is  a  metric  of  correction  stability  (lower 
is  better).  These  results  were  presented  at  the  2012  IEEE  Aerospace  conference. 
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Table  5.1:  Aggregate  AO  performance  gains  of  all  PCO-based  wavefront  reconstructors 
over  the  traditional  least  squares  reconstructor. 


Mean  Strehl 

Normalized  Strehl  Standard  Deviation 

0.04 

0.2% 

-2.4% 

0.5 

19.1% 

-50.6% 

1 

126.4% 

-57.6% 

•  Further  reconstructor  comparison  under  more  realistic  conditions.  It  turns  out 
LSPV+1  is  hard  to  beat  in  terms  of  performance  gain  relative  computations  required 
in  most  cases.  However,  as  the  Rytov  number  increases  above  roughly  0.5,  a  broader 
search  becomes  important  as  the  PDF  of  hopt  was  shown  to  broaden.  The  simplicity 
of  LSPV+4  with  very  similar  performance  to  the  more  advanced  reconstruction 
algorithms  proposed  by  Pellizzari  make  it  a  good  choice  in  high-Rytov  conditions. 

•  Identified  cause  of  fades  from  PCO-based  AO  and  demonstrated  a  potential  solution. 
When  in  closed-loop  operation  in  strong  turbulence,  PCO  algorithms  can  cause 
short-term  fades  as  hopt  varies.  These  drops  in  Strehl  are  not  caused  by  the  higher- 
order  control  law  as  much  as  different  tilt  components  associated  with  each  PCO 
rotational  realization.  The  tilt  is  removed  from  the  higher-order  corrections  and 
instantly  dumped  on  the  tracking  system;  causing  the  fades.  One  solution  was  found 
by  allowing  tilt  in  the  higher-order  corrections  with  conservative  controller  gains. 
This  prevented  the  fades  with  a  slight  tradeoff  in  mean  Strehl  ratio.  Unfortunately,  a 
realistic  system  with  slight  misalignments  would  require  a  zero-DC-gain  control  law 
to  implement  a  similar  solution.  Additionally,  this  puts  higher  stroke  requirements 
on  the  DM  when  tilt  is  being  included.  With  this  new  understanding  future  work  in 
control  design  could  help  mitigate  these  fades  based  on  the  application  requirements. 
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•  Demonstrated  and  developed  theory  on  improved  PCO  performance  in  weak 
turbulence,  low  SNR  regimes.  In  weak  turbulence,  where  no  branch  points  are 
present,  a  consistent  performance  improvement  of  the  PCO  reconstructors  over 
traditional  least  squares  was  found.  Although  the  difference  is  small,  this  was 
completely  unexpected  based  on  previous  concerns  of  noise  corrupting  the  PCO 
process.  A  new  theory  proposed  by  the  author  shows  in  high-noise  conditions, 
least  squares  reconstructors  can  introduce  phase  anomalies  to  the  AO  corrections. 
However,  the  PCO  process  removes  these  phase  discrepancies  introduced  by  the  least 
squares  reconstruction.  These  findings  will  be  presented  at  the  2012  SPIE  Optics  and 
Photonics  conference. 

•  Developed  a  set  of  theories  regarding  phase  anomalies  in  high  wavefront  gradient 
conditions.  When  localized  wavefront  phases  being  reconstructed  from  an  SRI  are 
greater  than  n  radians  between  adjacent  subapertures,  the  least  squares  reconstructor 
will  alias,  leading  to  a  phase  anomaly  in  the  reconstructed  field.  This  can  be  mitigated 
using  a  higher  density  of  SRI  subapertures,  at  the  cost  of  SNR.  Alternatively,  initial 
results  show  the  PCO  corrects  the  phase  anomalies  introduced  by  the  least  squares 
reconstruction. 

•  Formulated  an  optimized  data  flow  for  an  AO  system  using  PCO.  Started  by 
diagnosing  coupling  between  the  high-order  and  tracking  control  systems  when  using 
a  PCO-based  wavefront  reconstructor.  Since  PCO  reconstruction  is  not  a  single 
matrix  operation,  decoupling  tilt  from  the  higher-order  commands  required  a  new 
data  flow  process.  The  resulting  process,  illustrated  in  Fig.  3.4,  minimizes  matrix 
multiplies  while  maximizing  AO  effectiveness. 

•  Conducted  first  beam  projection  simulations  using  PCO  wavefront  reconstruction 
techniques.  Conducted  at  the  end  of  this  research,  this  establishes  a  foundation 
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for  future  efforts  to  investigate  these  novel  reconstruction  methods  with  respect 
to  beam  projection  applications.  Although  the  initial  simulation  results  identified 
needed  changes,  initial  observations  indicate  promise  over  traditional  least  squares 
reconstruction.  The  identified  flaws  have  since  been  fixed  and  follow-up  simulation 
results  will  be  presented  at  the  2012  DEPS  Beam  Control  Conference. 

•  Extensible  AFIT  AO  system  model  for  future  research.  The  model  is  configurable 
to  varying  levels  of  realism  and  AO  capabilities.  Allows  isolation  of  the  research 
problem  at  hand  without  the  difficulties  of  hardware  work  and  the  assumptions  of 
analytical  work.  This  includes  the  first  modeling  of  realistic  fiber  coupling  to  our 
specific  SRI  optical  setup  -  crucial  to  realistic  noise  modeling. 

5.3  Future  Work 

•  Hardware  verification  of  simulation  results.  Unfortunately,  the  DM  amplification 
circuit  board  failed  prior  to  beginning  hardware  testing  and  verification  of  the 
simulation  work  presented  here.  This  will  be  an  essential  step  forward  in  verifying 
all  theories  formulated  thus  far  regarding  PCO-based  wavefront  reconstruction  and 
control. 

•  Hardware  verifications  using  reed  atmospheric  turbulence.  Following  hardware 
verification  in  the  lab,  the  next  step  would  be  to  test  out  the  system  performance 
using  real  atmospheric  turbulence.  As  the  realism  is  increased,  new  challenges  are 
often  discovered  and  additional  validity  is  given  to  the  theory. 

•  Additional  verification  of  phase  anomaly  theories.  A  set  of  initial  simplified  tests 
were  made  to  form  the  theories  presented  on  phase  anomalies.  More  rigorous  testing 
to  relate  the  theories  to  more  realistic  scenarios  is  needed. 

•  Further  investigation  into  a  practical  tilt  coupling  between  the  control  systems  to 
reduce  fades.  A  more  in-depth  consideration  of  the  tilt  coupling  effects  between 
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the  tracking  system  and  a  higher-order  system  implementing  a  PCO  wavefront 
reconstructor  is  needed.  A  zero-DC-gain  controller  is  required  to  account  for  realistic 
component  misalignments,  which  was  not  implemented  here.  Another  alternative 
could  be  a  feed-forward  tilt-offloading  between  the  two  control  systems  as  h  is  varied. 

•  Investigate  adaptive  controlling  techniques  with  a  PCO-based  reconstructor.  Promis¬ 
ing  improvements  in  AO  performance  by  Gibson  [26]  and  Poyneer  [32]  using  adap¬ 
tive  AO  control  have  been  applied  towards  traditional  weak  turbulence  techniques. 
Investigation  as  to  whether  these  algorithms  would  work  with  the  PCO  non-unique 
rotational  phase  component  would  be  needed. 

•  Effects  of  SRI  performance  with  respect  to  platform  vibrations.  Moving  the  SRI  out 
of  a  lab  setup  could  reduce  the  fringe  stability  and  thus  the  fidelity  of  the  initial  SRI 
wrapped  measurements.  This  could  be  crudely  simulated  by  disturbing  the  lab  bench 
without  requiring  transportation  of  the  system  outside  the  lab. 

•  Further  previous  work  on  the  hybrid  SRI— Shack- Hartmann  WFS  using  the  simulation 
model  created  here.  The  simulation  model  developed  as  part  of  this  research  would 
be  an  excellent  and  easily  extended  platform  on  which  to  further  previous  work 
in  designing  an  effective  hybrid  WFS.  This  may  involve  dynamic  power  splitting, 
comparative  wavefront  reconstruction,  and  more  complex  control  techniques.  The 
range  of  investigations  and  potential  improvements  are  extensive,  and  this  would  be 
the  first  detailed  test  of  a  realistic  implementation  without  the  issues  involved  in  full 
hardware  integration. 

•  Further  thought  into  IWCL  and  CL  metrics  in  the  PCO  algorithms  that  operate 
over  the  residucd  phase  error,  as  opposed  to  the  fully  aberrated  field.  In  steady- 
state  closed  loop  operation,  the  wavefront  reconstructors  do  not  observe  the  fully 
aberrated  beacon  phase.  Rotational  phase  effects,  including  branch  cuts,  will  be  of 
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reduced  visibility  -  often  to  the  point  where  the  branch  cuts  are  below  the  threshold 
for  detecting  them  at  all.  Only  as  the  turbulence  evolves  in  time  are  features  visible, 
yet  still  not  completely  visible.  These  metrics  do  improve  PCO  performance  with 
respect  to  the  aliasing  phase  anomalies,  but  the  author  questions  whether  this  is  the 
best  way  to  evaluate  the  “goodness”  of  a  particular  reconstruction  when  operating  in 
closed  loop. 
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Appendix  A:  Detailed  Simulation  Results 


The  following  plots  show  the  compiled  simulation  results  for  each  reconstructor  over 
all  conditions  tested.  Each  grid-point  represents  a  simulated  one  second  AO  run.  Note  that 
the  simulations  with  higher  Greenwood  frequencies  will  run  over  many  more  Greenwood 
periods  than  the  lower  Greenwood  frequency  cases.  Unfortunately,  it  would  have  been 
computationally  prohibitive  to  run  all  of  the  low-Greenwood  cases  long  enough  to  hold  the 
number  of  Greenwood  periods  constant.  While  this  is  a  definite  compromise,  I  think  the 
usefulness  of  the  results  is  still  strong. 

The  left  column  is  the  mean  field-estimated  Strehl,  calculated  from  after  the  higher- 
order  control  system  is  enabled  through  the  end  of  each  simulation  run.  The  right  column 
shows  normalized  standard  deviation  of  the  prior  field-estimated  Strehl  ratio.  This  second 
statistic  provides  an  idea  for  the  system  stability,  where  lower  is  better.  The  color  axes  have 
been  fixed  for  each  section  below  for  easier  comparison. 
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Figure  A. 2:  LSPV+1  Reconstructor.  Weak  Turbulence:  cr2  =  0.04,  d / r0  =  1 
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Figure  A. 3:  LSPV+4  Reconstructor.  Weak  Turbulence:  cr2  =  0.04,  d/r0  =  1 
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Figure  A.4:  Histogram  Reconstructor.  Weak  Turbulence:  cr2  =  0.04,  d / r0  =  1 
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Figure  A.7:  Least-Squares  Reconstructor.  Medium  Turbulence:  cr~  =  0.5,  d/r0  =  1 
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Figure  A. 9:  LSPV+4  Reconstructor.  Medium  Turbulence:  crp  =  0.5,  cl / r()  =  1 


Mean  Strehl: 


Normalized  Standard  Deviation 


5  10  15 

SNR 


g  6°l 

3 


O 


10  15 

SNR 


0.4 

0.35 

0.3 


r 


Figure  A.  10:  Histogram  Reconstructor.  Medium  Turbulence:  <xp  =  0.5,  t//r0  =  1 
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Figure  A.  12:  Golden-Ratio  Reconstructor.  Medium  Turbulence:  crj  =  0.5,  d/r0=  1 
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Figure  A.  13:  Least-Squares  Reconstructor.  Strong  Turbulence:  crx  -  1,  t//r0  =  1 
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Figure  A.  15:  LSPV+4  Reconstructor.  Strong  Turbulence:  cr2x  =  l,  d/r0  =  1 
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Figure  A.  16:  Histogram  Reconstructor.  Strong  Turbulence:  cr2  =  1,  t//r0  =  1 
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Figure  A.  17:  Iterative  Probing  Reconstructor.  Strong  Turbulence:  cr^  =  1,  t//r0  =  1 
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Appendix  B:  PCO  Rotational  Phase  Tilt  Propagation 


This  appendix  presents  a  simple  simulation  tracing  a  rotational  wavefront  through  the 
higher-order  data  path.  It  is  related  to  how  the  tracking  system  perceives  the  wavefront  tilt 
both  before  and  after  the  higher-order  corrections  are  applied. 

The  simulation  starts  with  the  rotational  phase  shown  in  Fig.  B.l,  with  a  uniform 
irradiance  across  the  entire  field.  For  this  example,  the  wavefront  sensor  is  idealized, 
giving  no  sensing  errors  or  noise  effects.  The  least  squares  reconstruction  is  also  shown 
and  is  mostly  flat  as  discussed  in  Sec.  2.3.  Note  that  to  the  tracker,  the  rotational  phase 
does  not  have  any  tilt.  Figure  B.2  shows  the  tracker  spot  for  the  uncorrected  input  field.  It 
has  a  doughnut-like  shape,  but  it  is  centered  on  the  sensor  array. 

The  four  PCO  rotational  reconstructions  are  plotted  in  Fig.  B.3.  In  each  rotational 
realization,  a  tilt  component  is  defined  by  the  location  of  the  branch  cut.  It  is  important 
to  distinguish  this  tilt  component  from  the  tilt  component  viewable  by  the  tracker.  The  tilt 
component  with  respect  to  the  higher-order  system  is  defined  by  the  tilt-removal  operator 
from  Sec.  3.1.1.  To  this  operator,  each  rotational  realization  has  a  tilt  component  as  shown 
in  Fig.  B.4.  The  tilt-removal  operator  is  applied  to  all  four  rotational  realizations,  and  the 
resulting  outputs  are  applied  to  a  DM.  Sending  the  original  rotational  field  onto  each  DM 
realization,  and  then  looking  again  at  what  the  tracker  would  see  exposes  an  apparent  shift 
in  each  of  the  tracker  spots.  Figure  B.5  plots  each  tracker  output.  Note  that  the  spots  are 
no  longer  doughnut-like,  so  the  higher-order  phase  aberrations  were  correctly  compensated. 
However,  the  tilt  components  removed  from  each  rotational  realization,  not  only  are  visibly 
different,  but  they  clearly  propagate  to  the  tracking  system.  Unfortunately,  it  does  take  a 
finite  amount  of  time  for  the  tracking  system  to  react  to  these  changes.  That  is  the  cause  of 
the  fades  when  hopt  changes  in  a  high  Rytov  environment. 
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Figure  B.l:  A  hypothetical  rotational  phase  (left),  and  its  corresponding  least  squares 
reconstruction  (right). 


Tracker  Spot  for  Input  Beam  -  Centered,  Doughnut-like 


(Tpe) 


Figure  B.2:  The  tracker  image  of  the  input  field,  from  Fig.  B.l. 
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(frroti  h  —  7T  (frroti  h  37r/2 


Figure  B.3:  The  four  PCO  rotational  phase  realizations  from  an  LSPV+4  reconstructor 
over  the  rotational  field  from  Fig.  B.l. 
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Residual  Tilt,  h  —  0 
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Figure  B.4:  The  four  corresponding  tilt  components  from  each  rotational  phase  realization 
of  Fig.  B.3. 


117 


Tracker  Spot  for  h  =  0 


Tracker  Spot  for  h  —  it/ 2 


(Tpe) 


Tracker  Spot  for  h  —  it 


(Tpe) 


Tracker  Spot  for  h  =  37r/2 


1 

r  i14 

12 

10 
8 
6 
4 

2 
0 

(Tpe) 


Figure  B.5:  The  four  resulting  tracker  measurements  after  the  higher-order  corrections 
from  each  rotational  phase  realization  of  Fig.  B.3. 
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